perm filename V2D.IN[1,DEK] blob
sn#285515 filedate 1977-05-28 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00015 PAGES
C REC PAGE DESCRIPTION
C00001 00001
C00002 00002 folio 142 galley 1
C00014 00003 folio 144 galley 2
C00025 00004 folio 147 galley 3
C00041 00005 folio 151 galley 4
C00054 00006 folio 156 galley 5
C00069 00007 folio 161 galley 6
C00080 00008 folio 165 galley 7
C00093 00009 folio 173 galley 8
C00106 00010 folio 176 galley 9
C00120 00011 folio 178 galley 10
C00132 00012 folio 182 galley 11
C00145 00013 folio 186 galley 12
C00157 00014 folio 190 galley 13
C00173 00015 folio 194 galley 14
C00185 ENDMK
C⊗;
folio 142 galley 1
0 {U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addison-W
1 esley)|9f. 142|9ch. 3|9g. 1c|'{A6}!|9|7In principle,
7 any of these other random quantities may be obtained
16 from the uniform deviates |εU|β0, U|β1, U|β2,|4.|4.|4.|4.
23 |πPeople have devised a number of important ``random
31 tricks'' which may be used to perform these manipulations
40 e∃ciently on a computer, and a study of these
49 techniques also gives some insight into the proper
57 use of random numbers in any Monte Carlo application.|'
66 !|9|7It is conceivable that someday someone will
73 invent a random-number generator which produces
79 some of these other random quantities |εdirectly,
86 |πinstead of getting them indirectly via the
93 uniform distribution. Except for the ``random
99 bit'' generator described in Section 3.2.2, no
106 direct methods have so far proved to be practical.|'
115 !|9|7The discussion in the following section
121 assumes the existence of a random sequence of
129 uniformly distributed real numbers between zero
135 and one. A new uniform deviate is generated whenever
144 we need it. These numbers are normally represented
152 in a computer word with the decimal point assumed
161 at the left.|'{A18}|∨3|∨.|∨4|∨.|∨1|∨.|9|∨N|∨u|∨n|∨e|∨r|∨i|∨c
164 |∨a|∨l|4|4|∨D|∨i|∨s|∨t|∨r|∨i|∨b|∨u|∨t|∨i|∨o|∨n|∨s|'
165 {A6}This section summarizes the best techniques
171 known for producing numbers from various important
178 distributions. Many of the methods were originally
185 suggested by John von Neumann in the early 1950's,
194 and these have been gradually improved upon by
202 other people, notably George Marsaglia, J. H.
209 Ahrens, and U. Dieter.|'{A12}|≡A|≡.|9|≡R|≡a|≡n|≡d|≡o|≡m
214 |≡c|≡h|≡o|≡i|≡c|≡e|≡s |≡f|≡r|≡o|≡m |≡a |≡_|≡n|≡i|≡t|≡e
218 |≡s|≡e|≡t|≡.|9The simplest and most common type
224 of distribution required in practice is a random
232 |εinteger. |πAn integer between 0 and 7 can be
241 extracted from three bits of |εU |πon a binary
250 computer; in such a case, these bits should be
259 extracted from the |εmost signi⊂cant (|πleft-hand)
265 part of the computer word, since in many random-number
274 generators the least signi_cant bits are not
281 su∃ciently random. (See the discussion of this
288 in Section 3.2.1.1.)|'!|9|7In general, to get
295 a random integer |εX |πbetween 0 and |εk|4α_↓|41,
303 |πwe can |εmultiply |πby |εk, |πand let |εX|4α=↓|4|"lkU|"L.
311 |πOn |¬m|¬i|¬x, we would write|'{A9}|¬l|¬d|¬a!|¬u|;
317 |¬m|¬u|¬l!|¬k|;{A9}and after these two instructions
323 have been executed the desired integer will appear
331 in register A. If a random integer between 1
340 and |εk |πis desired, we add one to this result.
350 [The instruction ``|¬i|¬n|¬c|¬a |¬1'' would follow
356 (1).]|'!|9|7This method gives each integer with
363 nearly equal probability. (There is a slight
370 error because the computer word size is _nite;
378 see exercise 2. The error is quite negligible
386 if |εk |πis small, for example, if |εk/m|4|π|¬W|41/10000.)
394 In a more general situation we might want to
403 give di=erent weights to di=erent integers. Suppose
410 that the value |εX|4α=↓|4x|β1 |πis to be obtained
418 with probability |εp|β1, X|4α=↓|4x|β2 |πwith
423 probability |εp|β2,|4.|4.|4.|4, |πand |εX|4α=↓|4x|βk
427 |πwith probability |εp|βk. |πWe can generate
433 a uniform number |εU |πand let|'{A9}|h|ε|∂X|4α=↓|4|5|6|∂x|βk
439 ,!!|∂|πif!!|εp|β1|4α+↓|4p|β2|4α+↓|4.|4.|4.|4α+↓|4p|βk|βα_↓|β
439 1|4|¬Q|4U|4|¬Q|41.|E|n|;|L|Lx|β1,|L|πif!!|ε0|4|¬E|4|εU|4|¬W|
440 4p|β1;>{A2}{A2}|L{A6}|εX|4α=↓|4{B6}|5|6|Lx|β2,|L|πif!!|εp|β1
441 |4|¬E|4U|4|¬W|4p|β1|4α+↓|4p|β2;>|L|L|¬O|4|¬O|4|¬O>
443 |L|L|εx|βk,|L|πif!!|εp|β1|4α+↓|4p|β2|4α+↓|4|¬O|4|¬O|4|¬O|4α+
443 ↓|4p|βk|βα_↓|β1|4|¬E|4U|4|¬W|41.>{A9}|π(Note
445 that |εp|β1|4α+↓|4p|β2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βk|4α=↓|4
446 1.)|'!|9|7There is a ``best possible'' way to
454 do the comparisons of |εU |πagainst various values
462 of |εp|β1|4α+↓|4p|β2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βs,
464 |πas implied in (2); this situation is discussed
472 in Section 2.3.4.5. Special cases can be handled
480 by more e∃cient methods; for example, to obtain
488 one of the eleven values 2, 3,|4.|4.|4.|4, 12
496 with the respective probabilities |f1|d336|),
501 |f2|d336|),|4.|4.|4.|4, |f6|d336|),|4.|4.|4.|4,
503 |f2|d336|), |f1|d336|), we could compute two
509 independent random integers between 1 and 6 and
517 add them together.|'!|9|7However, none of the
524 above techniques is really the fastest way to
532 select |εx|β1,|4.|4.|4.|4,|4x|βk |πwith the correct
537 probabilities. A method which is considerably
543 more e∃cient in most cases, at a slight increase
552 in storage space, is explained below in connection
560 with the ``rectangle-wedge-tail'' method and
565 further illustrated in exercises 20 and 21. See
573 also exercise 28.|'{A12}|≡B|≡.|9|≡G|≡e|≡n|≡e|≡r|≡a|≡l
577 |≡m|≡e|≡t|≡h|≡o|≡d|≡s |≡f|≡o|≡r |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o|≡u|
579 ≡s |≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡i|≡o|≡n|≡s|≡.|9|4The
581 most general real-valued distribution may be
587 expressed in terms of its ``distribution function''
594 |εF(x), |πwhich speci_es the probability that
600 a random quantity |εX |πwill be less than or
609 equal to |εx;|'{A6}{A3}F(x)|4α=↓|4|πprobability(|εX|4|¬E|4x)
612 .|E|;(3)|?{A9}|πThis function always increases
618 monotonically from zero to one:|'{A9}|εF(x|β1)|4|¬E|4F(x|β2)
623 ,!!|πif!!|εx|β1|4|¬E|4x|β2;!!F(|→α_↓|¬X)|4α=↓|40,!!F(|→α+↓|¬
623 X)|4α=↓|41.|E|;(4)|?{A9}|πExamples of distribution
628 functions are given in Section 3.3.1, Fig. 3.
636 If |εF(x) |πis continuous and strictly increasing
643 (so that |εF(x|β1)|4|¬W|4F(x|β2) |πwhen |εx|β1|4|¬W|4x|β2),
648 |πit takes on all values between zero and one,
657 and there is an |εinverse function F|gα_↓|g1(y)
664 |πsuch that, for 0|4|¬W|4|εy|4|¬W|41,|'{A9}y|4α=↓|4F(x)!!|πi
668 f|6and|6only|6if!!|εx|4α=↓|4F|gα_↓|g1(y).|E|;
669 (5)|?{A9}|πA general way to compute a random
677 quantity |εX |πwith the continuous, strictly
683 increasing distribution |εF(x) |πis to set|'{A9}|εX|4α=↓|4F|
689 gα_↓|g1(U).|E|;(6)|?{A9}|πFor the probability
694 that |εX|4|¬E|4x |πis the probability that |εF|gα_↓|g1(U)|4|
700 ¬W|4x, |πi.e., the probability that |εU|4|¬W|4F(x),
706 |πand this is |εF(x).|'!|9|7|πThe problem now
713 reduces to one of numerical analysis, namely
720 to _nd good methods for evaluating |εF|gα_↓|g1(U)
727 |πto the desired accuracy. Numerical analysis
733 lies outsi{U0}{H10L12M29}W58320#Computer Programming|9(Knuth
folio 144 galley 2
735 /Addison-Wesley)|9f. 144|9ch. 3|9g. 2c|'{A6}!|9|7In
740 the _rst place, if |εX|β1 |πis a random variable
749 having the distribution |εF|β1(x) |πand if |εX|β2
756 |πis an independent random variable with the
763 distribution |εF|β2(x), |πthen|'{A9}|∂max(|εX|β1,|4X|β2)!!|∂
766 |πhas|6the|6distribution!!|εF|β1(x)F|β2(x),|h(x)|4α_↓|4F|β1(
766 x)F|β2(x).|n|;{A4}|π|Lmin(|εX|β1,|4X|β2)|L|πhas|6the|6distri
767 bution!!|εF|β1(x)|4α+↓|4F|β2(x)|4α_↓|4F|β1(x)F|β2(x).>
768 {B20}(7)|E|?{A20}{A9}|π(See exercise 4.) For
773 example, the uniform deviate |εU |πhas the distribution
781 |εF(x)|4α=↓|4x, |πfor 0|4|¬W|4x|4|¬E|41; |πif
785 |εU|β1, U|β2,|4.|4.|4.|4,|4U|βt |πare independent
789 uniform deviates, then max(|εU|β1, U|β2,|4.|4.|4.|4,|4U|βt)
794 |πhas the distribution function |εF(x)|4α=↓|4x|gt,
799 0|4|¬W|4x|4|¬E|41. |πThis is the basis of the
806 ``maximum of |εt'' |πtest given in Section 3.3.2.
814 Note that the inverse function in this case is
823 |εF|gα_↓|g1(y)|4α=↓|4{H11}|¬H{H10}|v4y|). |πIn
825 the special case |εtα=↓|π2, we see therefore
832 that the two formulas|'{A9}|εX|4α=↓|4{H11}|¬H{H10}U|)!!|πand
836 !!|εX|4α=↓|4|πmax(|εU|β1,|4U|β2)|E|;(8)|?{A9}|πwill
839 give equivalent distributions to the random variable
846 |εX, |πalthough this is not obvious at _rst glance.
855 We need not take the square root of a uniform
865 deviate.|'!|9|7The number of tricks like this
872 is endless: any algorithm which employs random
879 numbers as input will give a random quantity
887 with |εsome |πdistribution as output. The problem
894 is to _nd general methods for constructing the
902 algorithm, given the distribution function of
908 the output. Instead of discussing such methods
915 in purely abstract terms, we shall study how
923 they can be applied in important cases.|'{A12}|≡C|≡.|9|≡T|≡h
930 |≡e |≡n|≡o|≡r|≡m|≡a|≡l |≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡i|≡o|≡n|
932 ≡.|9|4Perhaps the most important nonuniform,
937 continuous distribution is the so-called |εnormal
943 distribution with mean zero and standard deviation
950 one|*/:|'|\{A9}F(x)|4α=↓|4|(1|d2{H11}|¬H{H10}|v22|≤p|)|)|4|↔j
951 |urx|)α_↓|¬X|)|4e|gα_↓|gt|i2|g/|g2|4dt.|E|;(9)|?
953 {A9}|πThe signi_cance of this distribution was
959 indicated in Section 1.2.10. Note that the inverse
967 function |εF|gα_↓|g1 |πis not especially easy
973 to compute; but we shall see that several other
982 techniques are available.|'{A12}(|*/|↔O|\)|9|εThe
986 polar method, |πdue to G. E. P. Box, M. E. Muller,
997 and G. Marsaglia.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
1001 |≡P|9(|εPolar method for normal deviates).|9|πThis
1006 algorithm calculates two independent normally
1011 distributed variables, |εX|β1 |πand |εX|β2.|'
1016 {A3}{I1.8H}|π|≡P|≡1|≡.|5|5[Get uniform variables.]
1019 Generate two independent random variables, |εU|β1,
1025 U|β2, |πuniformly distributed between zero and
1031 one. Set |εV|β1|5|¬L|52U|β1|4α_↓|41, V|β2|5|¬L|52U|β2|4α_↓|4
1034 1. (|πNow |εV|β1 |πand |εV|β2 |πare uniformly
1041 distributed between |→α_↓1 and |→α+↓1. On most
1048 computers it will be preferable to have |εV|β1
1056 |πand |εV|β2 |πrepresented in ⊗oatine-point form
1062 at this point.)|'{A3}|≡P|≡2|≡.|9[Compute |εS.]
1067 |πSet |εS|5|¬L|5V|=|β1|g2|4α+↓|4V|=|β2|g2.|'{A3}|≡P|≡3|≡.|9[
1069 |πIs |εS|4|¬R|41?] |πIf |εS|4|¬R|41, |πreturn
1074 to step P1. (Steps P1 through P3 are executed
1083 1.27 times on the average, with a standard deviation
1092 of 0.587; see exercise 6.)|'{A3}|≡P|≡4|≡.|9[Compute
1098 |εX|β1, X|β2.] |πSet |εX|β1, X|β2 |πaccording
1104 to the following two equations:|'{A9}!!|εX|β1|4α=↓|4V|β1|↔H|
1109 (|→α_↓2|4|πln|4|εS|d2S|),!!X|β2|4α=↓|4V|β2|↔H|(|→α_↓2|4|πln|
1109 4|εS|d2S|)|4.|E|;(10)|?{A9}{IC}|πThese are normally
1114 distributed variables desired.|'{A12}!|9|7To
1118 prove the validity of this method, we use elementary
1127 analytic geometry and calculus: If |εS|4|¬W|41
1133 |πin step P3, the point in the plane with Cartesian
1143 coordinates (|εV|β1,|4V|β2) |πis a |εrandom point
1149 uniformly distributed inside the unit circle.
1155 |πTransforming to polar coordinates |εV|β1|4α=↓|4R|4|πcos|4|
1159 ≤], |εV|β2|4α=↓|4R|4|πsin|4|≤], |πwe _nd |εS|4α=↓|4R|g2,
1164 X|β1|4α=↓|4{H11}|¬H{H10}|v4|→α_↓2|4|πln|4|εS|)|4|πcos|4|≤],
1165 |εX|β2|4α=↓|4{H11}|¬H{H10}|v4|→α_↓2|4|πln|4|εS|)|4|πsin|4|≤]
1165 . Using also the polar coordinates |εX|β1|4α=↓|4R|¬S|4|πcos|
1171 4|≤]|¬S, |εX|β2|4α=↓|4R|¬S|4|πsin|4|≤]|¬S, we
1174 _nd that |≤]|¬S|4α=↓|4|≤] and |εR|¬S|4α=↓|4{H11}|¬H{H10}|→α_
1178 ↓2|4|πln|4|εS|). |πIt is clear that |εR|¬S |πand
1185 |≤]|¬S are independent, since |εR |πand |≤] are
1193 independent inside the unit circle. Also, |≤]|¬S
1200 is uniformly distributed between 0 and |ε2|≤p;
1207 |πand the probability that |εR|¬S|4|¬E|4r |πis
1213 the probability that |→α_↓2|4|πln|4|εS|4|¬E|4r|g2,
1217 |πi.e., the probability that |εS|4|¬R|4e|gα_↓|gr|i2|g/|g2.
1222 |πThis equals 1|4α_↓|4|εe|gα_↓|gr|i2|g/|g2, |πsince
1226 |εS|4α=↓|4R|g2 |πis uniformly distributed between
1231 zero and one. The probability that |εR|¬S |πlies
1239 between |εr |πand |εr|4α+↓|4dr |πis therefore
1245 the derivative of |ε1|4α_↓|4e|gα_↓|gr|i2|g/|g2,
1249 |πnamely, |εre|gα_↓|gr|i2|g/|g2|4dr. |πSimilarly,
1252 the probability that |≤]|¬S lies between |ε|≤u
1259 |πand |ε|≤u|4α+↓|4d|≤u |πis (|ε1/2|≤p)|4d|≤u.
1263 |πTherefore the probability that |εX|β1|4|¬E|4x|β1
1268 |πand that |εX|β2|4|¬E|4x|β2 |πis|'{A9}{M30}|↔j|ur|)|¬T(|εr,
1272 |≤u)|3|¬G|3r|4|πcos|4|ε|≤u|¬Ex|β1,r|4|πsin|4|ε|≤u|¬Ex|β2|¬Y|
1272 )|4|(1|d22|≤p|)|4e|gα_↓|gr|i2|g/|g2|4r|4dr|4d|≤u|4α=↓|4|(1|d
1272 22|≤p|)|4|↔j|ur|)|¬T(x,y)|3|¬G|3x|¬Ex|β1,y|¬Ex|β2|¬Y|)|4e|gα
1272 _↓|g(|gx|i2|gα+↓|gy|i2|g)|g/|g2|4dx|4dy|'{A4}{M29}α=↓|4|↔a|↔
1273 K|(1|d22|≤p|)|4|↔j|urx|β1|)α_↓|¬X|)|4e|gα_↓|gx|i2|g/|g2|4dx|
1273 ↔s|↔a|↔K|(1|d22|≤p|)|4|↔j|urx|β2|)α_↓|¬X|)|4e|gα_↓|gy|i2|g/|
1273 g2|4dy|↔s.!(11)|?{A9}|πThis proves |εX|β1 |πand
1278 |εX|β2 |πare independent and normally distributed,
1284 as desired.|'{A12}|ε(|*/|↔P)|9|\Marsaglia's rectangle-wedge-t
1287 ail method.|9|4|πIn this method we use the distribution|'
1295 {A9}|εF(x)|4α=↓|4|↔K|(2|d2|≤p|)|4|↔j|urx|)0|)|4e|gα_↓|gt|i2|
1295 g/|g2|4dt!!x|4|¬R|40,|E|;(12)|?{A9}|πso that
1299 |εF(x) |πgives the distribution of the |εabsolute
1306 value |πof a normal deviate. After |εX |πhas
folio 147 galley 3
1314 {U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addison-W
1315 esley)|9f. 147|9ch. 3|9f. 3c|'{A6}!|9|7The rectangle-wedge-t
1320 ail approach is based on several important general
1328 techniques which we shall explore as we develop
1336 the algorithm. the _rst key idea is to regard
1345 |εF(x) |πas a |εmixture |πof several other functions,
1353 namely to write|'{A9}|εF(x)|4α=↓|4p|β1F|β1(x)|4α+↓|4p|β2F|β2
1356 (x)|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βnF|βn(x),|E|'
1357 (13)|?{A9}|πwhere |εF|β1, F|β2,|4.|4.|4.|4,|4F|βn
1361 |πare appropriate distributions and |εp|β1, p|β2,|4.|4.|4.|4
1366 ,|4p|βn |πare nonnegative probabilities which
1371 sum to 1. If we generate a random variable |εX
1381 |πby choosing distribution |εF|βj |πwith probability
1387 |εp|βj, |πit is easy to see that |εX |πwill have
1397 distribution |εF |πoverall. Some of the distributions
1404 |εF|βj(x) |πmay be rather di∃cult to handle,
1411 even harder than |εF |πitself, but we can usually
1420 arrange things so that the probability |εP|βj
1427 |πis very small in this case. Most of the distributions
1437 |εF|βj(x) |πwill be quite easy to accommodate,
1444 since they will be trivial modi_cations of the
1452 uniform distribution. The resulting method yields
1458 an extremely e∃cient program, since its |εaverage
1465 |πrunning time is very small.|'!|9|7It is easier
1473 to understad the method if we work with the |εderivatives
1483 |πof the distributions instead of the distributions
1490 themselves. Let|'{A6}{A3}|εf|1(x)|4α=↓|4F|¬S(x),!!f|βj(x)|4α
1492 =↓|4F|=|βj|¬S(x);|;{A9}|πthese are called the
1497 |εdensity |πfunctions of the probability distributions.
1503 Equation (13) becomes|'{A9}|εf(x)|4α=↓|4p|β1f|β1(x)|4α+↓|4p|
1506 β2f|β2(x)|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4p|βnf|βn(x).NSS|;
1507 *?|εf(x)|4α=↓|4p|β1f|β1(x)|4α+↓|4p|β2f|β2(x)|4α+↓|4|¬O|4|¬O|4
1507 |¬O|4α+↓|4p|βnf|βn(x).|E|;(14)|?{A9}|πEach |εf|βj(x)
1511 |πis |→|¬R0, and the total area under the graph
1520 of |εf|βj(x) |πso there is a convenient graphical
1528 way to display the relation (14): The area under
1537 |εf(x) |πis divided into |εn |πparts, with the
1545 part corresponding to |εf|βj(x) |πhaving area
1551 |εp|βj. |πSee Fig. 9, which illustrates the situation
1559 in the case of interest to us here, with |εf(x)|4α=↓|4F|¬S(x
1568 )|4α=↓|4{H11}|¬H{H10}|v4(2/|≤p)|)|εe|gα_↓|gx|i2|g/|g2;
1569 |πthe area under this curve has been divided
1577 into |εn|4α=↓|4|π37 parts. There are 12 rather
1584 large rectangels, which represent |εp|β1f|β1(x),|4.|4.|4.|4,
1588 p|β1|β2f|β1|β2(x); |πthere are 12 skinny little
1595 rectangles perched on top of these, which represent
1603 |εp|β1|β3f|β1|β3(x),|4.|4.|4.|4, p|β2|β4f|β2|β4(x);
1605 |πthere are 12 wedge-shaped pieces, which represent
1612 |εp|β2|β5f|β2|β5(x),|4.|4.|4.|4, p|β3|β6f|β3|β6(x);
1614 |πand the remaining part |εp|β3|β7f|β3|β7(x)
1619 |πis the ``tail,'' namely the entire graph of
1627 |εf(x) |πfor |εx|4|¬R|43.|'{A6}{H9L11}|π|≡F|≡i|≡g|≡.|7|≡9|≡.
1630 |9|4The density function divided into 37 parts.
1637 The area of each part represents the average
1645 number of times a random number with that density
1654 is to be computed.|;{A6}{H10L12}!|9|7The rectangular
1660 parts |εf|β1(x),|4.|4.|4.|4,|4f|β1|β2(x) |πrepresent
1663 |εuniform distributions. |πFor example, |εf|β3(x)
1668 |πrepresents a random variable uniformly distributed
1674 between |f14Fd32*?|) and |f3|d34|). We want to
1681 make it possible to determine the probabilities
1688 |εp|β1, p|β2,|4.|4.|4.|4, p|β1|β2 |πe∃ciently
1692 as the normal deviate is being generated, so
1700 (assuming we have a binary computer#the procedure
1707 for a decimal computer is similar) let us take
1716 each of these to be multiples of, say, |f1|d3256|).
1725 This means we want the area under |εp|βjf|βj(x)
1733 |πto be a multiple of |f1|d3256|), so the altitude
1742 of |εp|βjf|βj(x) |πis taken to be |"l64|εf(j/4)|"L/64,
1749 |πfor |ε1|4|¬E|4j|4|¬E|412. |πThe following table
1754 gives the resulting probabilities:|'{A9}|∂!!|∂!!|∂!!|∂!!|∂!!
1758 |∂!!|∂!!|∂!!|∂!!|∂!!|∂!!|∂!!|∂!!!!!!!!|9|∂!!|5|∂|E|;
1759 |>|εp|β1|;p|β2|;p|β3|;p|β4|;p|β5|;p|β6|;p|β7|;
1767 p|β8|;p|β9|;p|β1|β0|;p|β1|β1|;p|β1|β2|;(p|β1|β3|4α+↓|4|¬O|4|
1772 ¬O|4|¬O|4α+↓p|β3|β7)|;{A7}(15)|?{B7}>{A2}|>|f49|d3256|)|;
1777 |f45|d3256|)|;|f38|d3256|)|;|f30|d3256|)|;|f23|d3256|)|;
1781 |f16|d3256|)|;|f11|d3256|)|;|f6|d3256|)|;|f4|d3256|)|;
1785 |f2|d3256|)|;|f1|d3256|)|;|f0|d3256|)|;|f31|d3256|).|;
1789 >{A9}|πIt follows that the uniform distributions
1796 |εF|β1,|4.|4.|4.|4,|4F|β1|β2 |πmay be used |εp|β1|4α+↓|4|¬O|
1800 4|¬O|4|¬O|4α+↓ p|β1|β2|4|¬V|487.9 |πpercent of
1804 the time; about 88 percent of the area under
1813 |εf(x) |πis taken up by the large rectangles.|'
1821 !|9|7An important technique has been given by
1828 Marsaglia for e∃cient selection of the 13 alternatives
1836 represented in (15). Given a random binary integer
1844 between 0 and 255, whose binary representation
1851 is |εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8, |πwe procees
1855 as follows:|'{A9}{M29}|πif!|ε!|20|4|¬E|4b|β1b|β2b|β3b|β4|4|¬
1857 W|410,|'{A2}|πif!|ε|5|540|4|¬E|4|εb|β1b|β2b|β3b|β4b|β5b|β6|4
1858 |¬W|452,|'{A2}|πif!208|4|¬E|4|εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b
1859 |β8|4|¬W|4225,|'{A2}|πif!|ε225|4|¬E|4b|β1b|β2b|β3b|β4b|β5b|β
1860 6b|β7b|β8,|'{B54}{I13.9L}|πlet!|εX|4α=↓|4A[b|β1b|β2b|β3b|β4]
1861 |4α+↓|4|f1|d34|)U;|'{A2}|πlet!|εX|4α=↓|4B[b|β1b|β2b|β3b|β4b|
1862 β5b|β6]|4α+↓|4|f1|d34|)U;|'{A2}|πlet!|εX|4α=↓|4C[b|β1b|β2b|β
1863 3b|β4b|β5b|β6b|β7b|β8]|4α+↓|4|f1|d34|)U;|'{A2}|πgenerate
1865 |εX |πusing one of the distributions |εF|β1|β2,|4.|4.|4.|4,
1872 F|β3|β7, |πas described below.|'{A9}{IC}{M29}|πHere
1877 |εA, B, |πand |εC |πare auxiliary tables which
1885 appear in Table 1, and |εU |πis a uniform deviate
1895 between zero and one.|'{A9}{H9L11M29}|∂!!!!|6|6|∂!!!!!|∂!!!!
1899 !|9|∂!!!!!|9|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂|E|;{H8L10}|∨T|∨a|∨b|∨
1900 l|∨e|4|4|∨1|;{A6}{H9L11}|>|εA[0]|4α=↓|40|'|9A[5]|4α=↓|4|f1|d
1903 32|)|'|9B[40]|4α=↓|4|f1|d34|)|'|9B[46]|4α=↓|4|f3|d34|)|'
1906 |9C[208]|4α=↓|40|'|9C[214]|4α=↓|41|'|9C[220]|4α=↓|4|f7|d34|)
1908 |'>|>A[1]|4α=↓|40|'|9A[6]|4α=↓|4|f1|d32|)|'|9B[41]|4α=↓|4|f1
1913 |d34|)|'|9B[47]|4α=↓|41|'|9C[209]|4α=↓|4|f1|d34|)|'
1916 |9C[215]>|>A[2]|4α=↓|40|'|9A[7]|4α=↓|4|f3|d34|)|'
1920 |9B[42]|4α=↓|4|f1|d34|)|'|9B[48]|4α=↓|4|f3|d32|)|'
1922 |9C[210]|4α=↓|4|f1|d32|)|'|9C[216]|4α=↓|41|'|9C[222]|4α=↓|4|
1924 f9|d34|)|'>|>A[3]|4α=↓|4|f1|d34|)|'|9A[8]|4α=↓|41|'
1929 |9B[43]|4α=↓|4|f1|d32|)|'|9B[49]|4α=↓|4|f3|d32|)|'
1931 |9C[211]|4α=↓|4|f1|d32|)|'|9C[217]|4α=↓|4|f3|d32|)|'
1933 |9C[223]|4α=↓|4|f9|d34|)|'>|>A[4]|4α=↓|4|f1|d34|)|'
1937 |9A[9]|4α=↓|4|f5|d34|)|'|9B[44]|4α=↓|4|f3|d34|)|'
1939 |9B[50]|4α=↓|4|f7|d34|)|'|9C[212]|4α=↓|4|f3|d34|)|'
1941 |9C[218]|4α=↓|4|f3|d32|)|'|9C[224]|4α=↓|4|f5|d32|)|'
1943 >|>|;|;|9B[45]|4α=↓|4|f3|d34|)|'|9B[51]|4α=↓|42|'
1949 |9C[213]|4α=↓|4|f3|d34|)|'|9C[219]|4α=↓|4|f3|d32|)|'
1951 >{B2}|J#>{A9}{H10L12M29}|πHere |εA, B, |πand
1957 |εC |πare auxiliary tables which appear in Table
1965 1, and |εU |πis a uniform deviate between zero
1974 and one.|'!|9|7To see why this procedure works,
1982 consider for example the case |εj|4α=↓|44. F|β4(x)
1989 |πis the uniform distribution between |f3|d34|)
1995 and 1, and |εX|4α=↓|4|f3|d34|)|4α+↓|4|f1|d34|)U
1999 |πhas this distribution. The probability that
2005 this distribution is used in the above method
2013 is |f1|d316|) times the number of appearances
2020 of ``|f3|d34|)'' in the |εA |πtable, plus |f1|d364|)
2028 times the number of its appearances in |εB, |πplus
2037 |f1|d3256|) times the number in |εC; |πthis equals
2045 |f1|d316|)|4α+↓|4|f3|d364|)|4α+↓|4|f2|d3256|)|4α=↓|4|f30|d32
2045 56|)|4α=↓|4|εp|β4, |πso |εF|β4(x) |πis selected
2050 with the right probability. Of course, it would
2058 be even faster to use the following method:|'
2066 {A9}{M32}``if!!|ε0|4|¬E|4b|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8|4|
2066 ¬W|4225,!!|πlet!!|εX|4α=↓|4T[b|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β
2066 8]|4α+↓|4|f1|d34|)U''|'{A9}{M29}|πin place of
2070 the _rst three tests in (16). Here |εT |πwould
2079 be a table with 225 entries, 49 of which are
2089 ``0'', 45 of which are ``|f1|d34|)'', 38 are
2097 ``|f1|d32|)'', etc. By comparison, the previous
2103 method uses only 39 table entries; and it is
2112 only a little slower, since the test ``|εb|β1b|β2b|β3b|β4b|β
2119 5b|β6|4|¬W|452?'' |πneeds to be made only three-eighths
2126 of the time, and the test ``|εb|β1b|β2b|β3b|β4b|β5b|β6b|β6b|
2132 β8|4|¬W|4225?'' needs to be made only three-sixteenths
2139 of the time. The savings of storage space in
2148 a method like (16) is even more signi_cant in
2157 the more general situation when there are more
2165 than 256 possibilities.|'{A6}|π|≡F|≡i|≡g|≡.|7|≡1|≡0|≡.|9|4Fr
2168 equency functions for which Algorithm L may be
2176 used to generate random numbers.|;{A6}{H10L12}!|9|7|πThe
2182 reader should pause at this point to study the
2191 above method carefully, until he understands
2197 why method (16) selects the distributions |εF|βj
2204 |πwith the probability |εp|βj, |πfor 1|4|¬E|4j|4|¬E|412.|'
2210 !|9|7The distributions |εF|β1|β3(x),|4.|4.|4.|4,|4F|β2|β4(x)
2212 |πare also uniform distributions, and they may
2220 be handled similarly. They have quite small probability,
2228 however (some number between 0 and |f1|d3256|)-,
2235 so it is most e∃cient to test for these probabilities
2245 in a part of our normal-deviate-generating alg{U0}{H10L12M29
folio 151 galley 4
2251 }w58320#Computer Programming|9(Knuth/Addison-Wesley)|9f.
2253 151|9ch. 3|9g. 4c|'{A6}!|9|7Now we turn to the
2261 question of computing random variables from the
2268 wedge-shaped distributions |εF|β2|β5(x),|4.|4.|4.|4,
2271 F|β3|β6(x). |πTypical cases are shown in Fig.
2278 10; when |εx|4|¬W|41, |πthe curved part is concabe
2286 downward, and when |εx|4|¬Q|41 |πit is concave
2293 upward, but in each case the curved part is reasonably
2303 close to a straight line, and it can be enclosed
2313 in two parallel lines as shown.|'{A12}!|9|7In
2320 order to handle these wedge-shaped distributions,
2326 we will rely on yet another general technique,
2334 von Neumann's so-called |εrejection method |πfor
2340 obtaining a complicated density from another
2346 one which ``encloses'' it. The polar method described
2354 above is a simple example of such an approach:
2363 Steps P1<P3 obtain a random point inside the
2371 unit circle by _rst generating a random point
2379 inside the unit circle by _rst generating a random
2388 point in a larger square, rejecting it and starting
2397 over again if the point was outside the circle.|'
2406 !|9|7The general rejection method is even more
2413 powerful than this. To generate a random variable
2421 |εX |πwith density |εf, |πlet |εg |πbe another
2429 probability density function such that |εf|1(t)|4|¬E|4cg(t)
2435 |πfor all |εt, |πwhere |εc |πis a constant. Now
2444 generate |εX |πaccording to density |εg, |πand
2451 also generate an independent uniform deviate
2457 |εU. |πIf |εU|4|¬R|4f|1(X)/cg(X), |πreject |εX
2462 |πand start again with another |εX |πand |εU.
2470 |πWhen the condition |εU|4|¬W|4f|1(X)/cg(X) |π_nally
2475 occurs, the resulting |εX |πwill have density
2482 |εf |πas desired. [Prgoof*?e resulting |εX |πwill
2489 have density |εf |πas desired. [Proof: |εX|4|¬E|4x
2496 |πwill occur with probability |εp(x)|4α=↓|4|↔j|urx|)0|)|4{H1
2500 2}({H10}g(t)|4dt|4|¬O|4f|1(t)/cg(t)}h12}){H10}|4α+↓|4gp(x),
2501 |πwhere |εq|4α=↓|4|↔j|ur1|)0|)|4{H12}({H10}g(t)|4dt|4|¬O|4(1
2502 |4α_↓|4f|1(t)/cg(t){H12}){H10}|4α=↓|41|4α_↓|41/c
2503 |πis the probability of rejection; hence |εp(x)|4α=↓|4|↔j|ur
2509 x|)0|)|4f(t)|4dt.]|'!|9|7|πThe rejection technique
2513 is most e∃cient when |εc |πis small, since there
2522 will be |εc |πiterations on the average before
2530 a value is accepted. (See exercise 6.) In some
2539 cases |εf|1(x)/cg(x) |πis hard to compute, we
2546 may know some bounding functions |εr(x)|4|¬E|4f|1(x)/cg(x)|4
2551 |¬E|4s(x) |πthat are much simpler, and the exact
2559 value of |εf(x)/cg(x) |πneed not be calculated
2566 unless |εr(x)|4|¬E|4U|4|¬W|4s(x). |πThe following
2570 algorithm solves the ``wedge'' problem by developing
2577 the rejection method still further.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i
2582 |≡t|≡h|≡m |≡L (|εNearly linear densities).|9|4|πThis
2587 algorithm may be used to generate a random variable
2596 |εX |πfor any distribution whose density |ε(x)
2603 |πsatis_es the following conditions (cf. Fig.
2609 10):|'{A9}{A8}(18)|E|?{B8}|εf|1(x)|4α=↓|40,!!|πfor!!|εx|4|¬W
2611 |4s!!|πand|6for!!|εx|4|¬Q|4s|4α+↓|4h;|;{A4}a|4α_↓|4b(x|4α_↓|
2612 4s)/h|4|¬E|4f|1(x)|4|¬E|4b|4α_↓|4b(x|4α_↓|4s)/h,!!|πfor!!|εs
2612 |4|¬E|4x|4|¬E|4s|4α+↓|4h.|;{A9}{I1.8H}|π|≡L|≡1|≡.|9[Get
2614 |εU|4|¬E|4V.] |πGenerate two independent random
2619 variables |εU, V, |πuniformly distributed between
2625 zero and one. If |εU|4|¬Q|4V, |πexchenge |εU|6|≠l|6V.|'
2632 {A3}|π|≡L|≡2|≡.|9[|πEasy case?] If |εV|4|¬E|4a/b,
2636 |πgo to L4.|'{A3}|≡L|≡3|≡.|9[Try again?] If |εV|4|¬Q|4U|4α+↓
2642 |4(1/b)|1f|1(s|4α+↓|4hU), |πgo back to step L1.
2648 (If |εa/b |πis close to 1, this step of the algorithm
2659 will not be necessary very often.)|'{A3}|≡L|≡4|≡.|9[Compute
2666 |εX.] |πSet |εX|5|¬L|5s|4α+↓|4hU.|'{IC}{A9}{H9L11}|π|≡F|≡i|≡
2669 g|≡.|7|≡1|≡1|≡.|9|4Region of ``acceptance'' in
2673 Algorithm L.|;{A9}{H10L12}!|9|7To prove that
2678 this algorithm is valid, we observe that when
2686 step L4 is reached, the point (|εU,|4V) |πis
2694 a random point in the area shaded in Fig. 11,
2704 namely, 0|4|¬E|ε|4U|4|¬E|4V|4|¬E|4U|4α+↓|4(1/b)|1f|1(s|4α+↓|
2705 4hU). |πConditions (18) ensure that|'{A9}|ε|(a|d2b|)|4|¬E|4U
2710 |4α+↓|4|(1|d2b|)|4f|1(s|4α+↓|4hU)|4|¬E|41.|;{A9}|πNow
2712 the probability that |εX|4|¬E|4s|4α+↓|4hx, |πfor
2717 0|4|¬E|4x|4|¬E|41, |πis the ratio of area to
2724 the left of the vertical line |εU|4α=↓|4x |πin
2732 Fig. 11 to the total area, namely,|'{A9}|↔j|ur|εx|)0|)|4|(1|
2739 d2b|)|4f|1(s|4α+↓|4hu)|4du|↔h|↔j|ur1|)0|)|4|(1|d2b|)|4f|1(s|
2739 4α+↓|4hu)|4du|4α=↓|4|↔j|ursα+↓hx|)s|)|4f|1(v)|4dv;|;
2740 {A9}|πtherefore |εX |πhas the correct distribution.|'
2746 {A9}{H9L11}|≡F|≡i|≡g|≡.|7|≡1|≡2|≡.|9|4The ``rectangle-wedge'
2747 ' algorithm for generating normal deviates.|;
2753 {A9}{H10L12}!|9|7With appropriate constants |εa|βj,
2757 b|βj, s|βj, |πAlgorithm L will take care of the
2766 wedge-shaped densities |εf|βj|βα+↓|β2|β4 |πof
2770 Fig. 9, for |ε1|4|¬E|4j|4|¬E|412. |πThe _nal
2776 distribution, |εF|β3|β7, |πneeds to be treated
2782 only about one time in 370; it is used whenever
2792 a result |εX|4|¬R|43 |πis to be given. Exercise
2800 11 shows that a standard rejection scheme can
2808 be used for this ``tail''; hence we are ready
2817 to consider the rectangle-wedge-tail procedure
2822 in its entirety:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
2826 |≡M (|εMarsaglia-MacLaren method for normal deviates).|9|4|π
2831 This algorithm uses several auxiliary tables,
2837 constructed as explained in the text, and examples
2845 appear in Tables 1 and 2. We assume that a binary
2856 computer is being used; the procedure for a decimal
2865 computer is similar.|'{A3}{H10L12}{I2.3H}|5|5|≡M|≡1|≡.|9[Get
2868 |εU.] |πGenerate a uniform random number |εU|4α=↓|4.b|β0b|β
2875 1b|β2|4.|4.|4.|4b|βm. (|πHere the |εb'|πs are
2880 the bits in the binary representation of |εU.
2888 |πFor reasonable accuracy, |εm |πshould be at
2895 least 24.) Set |ε|≤c|5|¬L|5b|β0. (|πLater, |ε|≤c
2901 |πwill be used to determine the sign of the result.)|'
2911 {A3}|5|5|≡M|≡2|≡.|9[Large rectangle?] If |εb|β1b|β2b|β3b|β4|
2914 4|¬W|410, |πwhere ``|εb|β1b|β2b|β3b|β4'' |πdenotes
2918 the binary integer |ε8b|β1|4α+↓|44b|β2|4α+↓|42b|β3|4α+↓|4b|β
2921 4, |πset |εX|5|¬L|5A[b|β1b|β2b|β3b|β4]|4α+↓|4.00b|β5b|β6|4.|
2923 4.|4.|4b|βm |πand go to M10. Otherwise, if |εb|β1b|β2b|β3b|β
2930 4b|β5b|β6|4|¬W|452, set |εX|5|¬L|5B[b|β1b|β2b|β3b|β4b|β5b|β6
2932 ]|4α+↓|4.00b|β7b|β8|4.|4.|4.|4b|βm |πand go to
2936 M10. Otherwise, if |εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8|4|¬W|4
2939 225, |πset |εX|5|¬L|5C[b|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8]|4α+
2941 ↓|4.00b|β9b|β1|β0|4.|4.|4.|4b|βm |πand go to
2945 M10.|'{A9}{IC}{H8L10}|∨T|∨a|∨b|∨l|∨e|4|4|∨2|;
2947 {A6}{H9L11}EXAMPLE|4|4OF|4|4TABLES|4|4USED|4|4WITH|4|4ALGORI
2947 THM|4|4Mα/↓|;{M27}{B2}|J#>|∂!!!!!|∂!!!!!|∂!!!!!!!|∂!!!!!!|9|
2949 ∂!!!!!!|∂!!!!!!|9|∂|E|;|>|εj|;S[j]|;P[j]|;Q[j]|;
2955 D[j]|;E[j]|;>{B2}|J#>|>|5|51|;0|;0.885|;0.881|;
2964 0.51|;16!|7|;>|>|5|52|;|f1|d34|)|;0.895|;0.885|;
2972 0.79|;|5|58!|7|;>|>|5|53|;|f1|d32|)|;0.910|;0.897|;
2980 0.90|;|5|55.33|;>|>|5|54|;|f3|d34|)|;0.929|;0.914|;
2988 0.98|;|5|54!|7|;>|>|5|55|;1|;0.945|;0.930|;0.99|;
2997 |5|53.08|;>|>|5|56|;|f5|d34|)|;0.960|;0.947|;
3004 0.99|;|5|52.44|;>|>|5|57|;|f3|d32|)|;0.971|;0.960|;
3012 0.98|;|5|52.00|;>|>|5|58|;|f7|d34|)|;0.982|;0.974|;
3020 0.96|;|5|51.67|;>|>|5|59|;2|;0.987|;0.982|;0.95|;
3029 1.43|;>|>10|;|f9|d34|)|;0.991|;0.989|;0.93|;|5|51.23|;
3038 >|>11|;|f5|d32|)|;0.994|;0.992|;0.94|;|5|51.08|;
3046 >|>12|;|f11|d34|)|;0.997|;0.996|;0.94|;|5|50.95|;
3054 >|>13|;3|;1.000|;>{B2}|J#>|Ha*?*?*?{U0}{H10L12M29}{I2.3H}W58320
folio 156 galley 5
3061 #Computer Programming|9(Knuth/Addison-Wesley)|9f.|4156|9ch.
3063 3|9g. 5c|'{A6}|5|5|≡M|≡3|≡.|9[Wedge or tail?]
3068 Find the |εsmallest |πvalue of |εj, 1|4|¬E|4j|4|¬E|413,
3075 |πfor which |ε.b|β1b|β2|4.|4.|4.|4b|βm|4|¬W|4P[j].
3078 |πIf |εj|4α=↓|413, |πgo to step M8.|'{A3}|5|5|≡M|≡4|≡.|9[Ski
3084 nny rectangle?] If |ε.b|β1b|β2|4.|4.|4.|4b|βm|4|¬W|4Q[j],
3088 |πgenerate a new unidorm random number |εU, |πset
3096 |εX|5|¬L|5S[j]|4α+↓|4|f1|d34|)U, |πand go to
3100 M10.|'{A3}|5|5|≡M|≡5|≡.|9[Get |εU|4|¬E|4V.] |πGenerate
3104 two new uniform deviates, |εU, V; |πif |εU|4|¬Q|4V,
3112 |πexchange |εU|5|≠m|5V. (|πWe are now performing
3118 Algorithm L.) Set |εX|5|¬L|5S[j]|4α+↓|4|f1|d34|)U.|'
3123 {A3}|5|5|π|≡M|≡6|≡.|9[|πEasy case?] If |εV|4|¬E|4D[j],
3127 |πgo to M10.|'{A3}|5|5|≡M|≡7|≡.|9[Another try?]
3132 If |εV|4|¬Q|4U|4α+↓|4E[j] (e|gα_↓|g(|gX|i2|gα_↓|gS|g[|gj|gα+
3134 ↓|g1|g]|i2|g)|g/|g2|4α_↓|41), |πgo back to step
3139 M5; otherwise go to M10. (This step is executed
3148 with low probability.)|'{A3}|5|5|≡M|≡8|≡.|9[Get
3152 |εU|g2|4α+↓|4V|g2|4|¬W|41.] |πGenerate two new
3156 independent uniform deviates, |εU, V; |πset |εW|5|¬L|5U|g2|4
3162 α+↓|4V|g2. |πIf |εW|4|¬R|41, |πrepeat this step.|'
3168 {A3}|5|5|≡M|≡9|≡.|9[Compute |εX|4|¬R|43.] |πSet
3171 |εT|4|¬L|4{H11}|¬H{H10}|v4(9|4α_↓|42|4|πln|4|εW)/W|).
3172 |πSet |εX|4|¬L|4U|4α⊗↓|4T. |πIf |εX|4|¬R|43,
3176 |πgo to M10. Otherwise set |εX|4|¬L|4V|4α⊗↓|4T.
3182 |πIf |εX|4|¬R|43, |πgo to M10. Otherwise go back
3190 to step M8. (The latter will occur about half
3199 the time this step is encountered.)|'{A3}|5|5|≡M|≡8|≡.|9[Get
3205 supertail deviate.] Generate two new independent
3212 uniform deviates, |εU, V, |πand set |εX|4|¬L|4{H11}|¬H{H10}|
3218 v49|4α_↓|42|4|πln|4|εV|).|'{A3}|5|5|π|≡M|≡9|≡.|9[|πReject
3220 ?] If |εUX|4|¬R|43, |πgo back to stop M8. (This
3229 will occur only about one-twelfth of the time.)|'
3237 {A3}|5|5|≡M|≡1|≡0|≡.|9[Attach sign.] If |ε|≤c|4α=↓|41,
3241 |πset |εX|4|¬L|4|→α_↓X.|'{IC}{A12}|π!|9|7This
3244 algorithm is a very pretty example of mathematical
3252 theory intimately interwoven with programming
3257 ingenuity#it is a _ne illustration of the art
3265 of computer programming.|'!|9|7Tables A, B, and
3272 C have already been described; the other tables
3280 required by Algorithm M are constructed as follows
3288 for 1|4|¬E|4|εj|4|¬E|412, |πin terms of the quantities
3295 |εa|βj, b|βj, |πand |εp|βj|βα+↓|β2|β4 |πde_ned
3300 in exercise 10:|'{A6}|εP[j]|4α=↓|4p|β1|4α+↓|4p|β2|4α+↓|4|¬O|
3303 4|¬O|4|¬O|4α+↓|4p|β1|β2|4α+↓|4(p|β1|β3|4α+↓|4p|β2|β5)|4α+↓|4
3303 |¬O|4|¬O|4|¬O|4α+↓|4(p|βj|βα+↓|β1|β2|4α+↓|4p|βj|βα+↓|β2|β4);
3303 |'P[13]|4α=↓|41; S[j]|4α=↓|4(j|4α_↓|41)/4; Q[j]|4α=↓|4P[j]|4
3306 α_↓|4p|βj|βα+↓|β2|β4;|'D[j]|4α=↓|4a|βj/b|βj;
3308 E[j]|4α=↓|4{H11}|¬H{H10}2/|≤p|)|4|πexp{H12}({H10}|→α_↓(j/4)|
3308 g2/2)/b|βjp|βj|βα+↓|β2|β4.|E|'(19)|?|πTable 2
3312 shows the values for our example to only a few
3322 signi_cant _gures, but in a computer program
3329 these quantities would all appear with full-word
3336 precision. Algorithm M requires a total of 101
3344 auxiliary table entries in all.|'!|9|7This method
3351 is extremely fast, since 88 percent of the time
3360 only the steps M1, M2, and M10 will need to be
3371 executed, and the other steps aren't terribly
3378 slow either. We have broken the range from 0
3387 to 3 into 12 parts in Fig. 9; if it were broken
3399 into more parts, say 48 parts, more table entries
3408 would be required, but the computation would
3415 stay in steps M1, M2, M10 over 97 percent of
3425 the time*3 Complete tables for both binary and
3433 decimal computers are given in the article by
3441 Marsaglia, MacLaren, and Bray, |εCACM|67 (1964),
3447 4<10; |πan additional trick, overlapping parts
3453 of Tables A, B, C and S, is employed there in
3464 order to save storage space. Further re_nements
3471 have been developed by Marsaglia, Anathanarayanax,
3477 and Paul, |εInf. Proc. Letters |≡4 (1976), |πto
3485 appear.|'{A12}|ε|*/(|↔L) |\Forsythe's method.|9|4|πAn
3489 amazingly simple technique for generating random
3495 deviates with a density of the general exponential
3503 form|'{A9}|εf|1(x)|4α=↓|4Ce|gα_↓|gh|g(|gx|g),!!|πfor!!|εa|4|
3504 ¬E|4x|4|¬W|4b;|6f|1(x)|4α=↓|40|6|πotherwise;|E|;
3505 (20)|?{A4}|ε0|4|¬E|4h(x)|4|¬E|41!!|πfor!!|εa|4|¬E|4x|4|¬W|4b
3506 ;|E|;(21)|?{A9}|πwas discovered by John von Neumann
3514 and G. E. Forsythe about 1950. The idea is based
3524 on the rejection method described earlier, letting
3531 |εg(x) |πbe the uniform distribution on [|εa,b):
3538 |πWe set |εX|4|¬L|4a|4α+↓|4(b|4α_↓|4a)U, |πwhere
3542 |εU |πis a uniform deviate, and then we want
3551 to accept |εX |πwith probability |εe|gα_↓|gh|g(|gX|g).
3557 |πThe latter operation could be done by testing
3565 |εe|gα_↓|gh|g(|gX|g) |πvs. |εV, |πor |εh(X) |πvs.
3571 |ε|→α_↓|πln |εV, |πwhen |εV |πis another uniform
3578 deviate, but the job can be done without applying
3587 any transcendental functions in the following
3593 interesting way. Set |εV|β0|4|¬L|4h(X), |πthen
3598 generate uniform deviates |εV|β1, V|β2,|4.|4.|4.
3603 |πuntil _nding some |εK|4|¬R|41 |πwith |εV|βK|βα_↓|β1|4|¬W|4
3608 V|βK. |πFor _xed |εV |πand |εk, |πthe probability
3616 that |εH(X)|4|¬R|4V|β1|4|¬R|4|¬O|4|¬O|4|¬O|4|¬R|4V|βk
3618 |πis |ε1/k*3 |πtimes the probability that max|4(|εV|β1,|4.|4.
3624 |4.|4,|4V|βk){H10}|4|¬E|4|εh(X), |πnamely |εh(X)|gk/k*3;
3627 |πhence the probability that |εK|4α=↓|4k |πis
3633 |εh(X)|gk|gα_↓|g1/(k|4α_↓|41)*3|4α_↓|4h(X)|gk/k*3,
3634 |πand the probability that |εK |πis odd is|'{A9}|ε|↔k|uc|)k|
3642 4|πodd|d|εk|¬R1|)|4|↔a|(h(X)|gk|gα_↓|g1|d2(k|4α_↓|41)*3|)|4α_
3642 ↓|4|(h(X)|gk|d2k*3|)|↔s|4α=↓|4e|gα_↓|gh|g(|gX|g).|;
3643 {A9}|πTherefore we reject |εX |πand try again
3650 if |εK |πis even; we accept |εX |πas a random
3660 variable with density (20) if |εK |πis odd. Note
3669 that we usually won't have to generate many |εV'|πs
3678 in order to determine |εK, |πsince the average
3686 value of |εK (|πgiven |εX) |πis |¬K|ur|)|εk|¬R0|)|4|πPr(|εK|
3692 4|¬Q|4k)|4α=↓|4|¬K|ur|)k|¬R0|)|4h(X)|gk/k*3|4α=↓|4e|gh|g(|gX|
3692 g)|4|¬E|4e.|'!|9|7|πForsythe realized some years
3697 later that this approach leads to an e∃cient
3705 method for calculating normal deviates, without
3711 the need for any |πauxiliary routines to calculate
3719 square roots or logarithms as in Algorithms P
3727 and M. His procedure, with an improved choice
3735 of intervals [|εa,b) |πdue to J. H. Ahrens and
3744 U. Dieter, can be summarized as follows.|'{A12}|≡A|≡l|≡g|≡o|
3751 ≡r|≡i|≡t|≡h|≡m |≡F (|εForsythe method for normal
3757 deviates).|9|4|πThis algorithm generates normal
3761 deviates on a binary computer, assuming approximately
3768 |εm|4α+↓|41 |πbits of accuracy. |πExtremely large
3774 values of |εX |πwhose probability is less than
3782 2|ε|gα_↓|gm|gα_↓|g1 |πwill never occur. The algorithm
3788 requires a table of values |εd|βj|4α=↓|4a|βj|4α_↓|4a|βj|4α_↓
3793 |41, |πfor |ε1|4|¬E|4j|4|¬E|4m|4α+↓|41, |πwhere
3797 |εa|βj |πis de_ned by the relation|'{A9}|ε|↔H|(2|d2|≤p|)|4|↔
3803 j|ur|¬X|)a|βj|)|4e|gα_↓|gt|i2|g/|g2|4dt|4α=↓|4|(1|d22|gj|)|4
3803 .|E|;(23)|?{A9}{I1.8H}|π|≡F|≡1|≡.|9[Get |εU.]
3807 |πGenerate a uniform random number |εU|4α=↓|4.b|β0b|β1|4.|4.
3812 |4.|4b|βm, |πwhere |εb|β0,|4.|4.|4.|4,|4b|βm
3815 |πdenote the bits in binary notation. Set |ε|≤c|4|¬L|4b|β0,
3823 j|4|¬L|41, |πand |εa|4|¬L|40.|'{A3}|π|≡F|≡2|≡.|9[Find
3827 _rst zero |εb|βj.]|9|πIf |εb|βj|4α=↓|41, |πset
3832 |εa|4|¬L|4a|4α+↓|4d|βj, j|4|¬L|4j|4α+↓|41, |πand
3835 |εrepeat |πthis step. (If |εj|4α=↓|4m|4α+↓|41,
3840 |πtreat |εb|βj |πas zero.)|'{A3}|≡F|≡3|≡.|9[Generate
3845 candidate.]|9(Now |εa|4α=↓|4a|βj|βα_↓|β1, |πand
3848 the current value of |εj |πoccurs with probability
3856 |→|¬V1/2|ε|gj. |πWe will generate |εX |πin the
3863 range [|εa|βj|βα_↓|β1,|4a|βj), |πwing the rejection
3868 method described above, with |εh(x)|4α=↓|4x|g2/2|4α=↓|4y|g2/
3872 2|4α+↓|4ay |πwhere |εy|4α=↓|4x|4α_↓|4a. |πExercise
3876 12 proves that |εh(x)|4|¬E|41 |πas required in
3883 (21).) Set |εY|4|¬L|4d|βj |πtimes |ε.b|βj|βα+↓|β1|4.|4.|4.|4
3887 b|βm |πand |εV|4|¬L|4(|f1|d32|)Y|4α+↓|4a)Y. (|πSince
3891 the average value of |εj |πis 2, there will usually
3901 be enough signi_cant bits in |ε.b|βj|βα+↓|β1|4.|4.|4.|4b|βm
3907 |πto provide decent accuracy. The calculations
3913 are readily done in _xed point arithmetic.)|'
3920 {A3}|≡F|≡4|≡. [Reject*3]|9|4Generate a uniform
3924 deviate |εU. |πIf |εV|4|¬W|4U, |πgo on to step
3932 F5. Otherwise set |εV |πto a new uniform deviate;
3941 and if now |εU|4|¬W|4V (|πi.e., if |εK|πis even
3949 in the discussion above), go back to F3, otherwise
3958 repeat s{U0}{H10L12M29}|πW58320#Computer Programming|9(Knuth
folio 161 galley 6
3960 /Addison-Wesley)|9f. 161|9ch. 3|9g. 6c|'{A6}!|9|7Values
3965 of |εd|βj |πfor |ε1|4|¬E|4j|4|¬E|447 |πappear
3970 in a paper by Ahrens and Dieter, |εMath. comp.
3979 |π|≡2|≡7 (1973), 927<937. This paper discusses
3985 re_nements of the algorithm which improve its
3992 speed at the expense of more tables. Algorithm
4000 F is attractive since it is almost as fast as
4010 Algorithm M and it is much easier to implement.
4019 The average number of uniform deviates required
4026 per normal deviate is 2.53947; Brent [(|εACM
4033 |π|≡1|≡7 (1974), 704<705] has shown how to reduce
4041 this number to 1.37446 at the expense of two
4050 subtractions and one division per uniform deviate
4057 saved.|'{A12}|ε|*/(|↔C)|9|\Variations of the normal
4062 distribution.|9|4|πSo far we have considered
4067 the normal distribution with mean zero and standard
4075 deviation one. If |εX |πhas this distribution,
4082 then|'{A9}|εY|4α=↓|4|≤m|4α+↓|4|≤sX|E|;(24)|?{A9}|πhas
4086 the normal distribution with mean |ε|≤m |πand
4093 standard deviation |ε|≤s. |πFurthermore, if |εX|β1
4099 |πand |εX|β2 |πare independent normal deviates
4105 with mean zero and standard deviation one, and
4113 if|'{A9}|εY|β1|4α=↓|4|≤m|β1|4α+↓|4|≤s|β1X|β1,!!Y|β2|4α=↓|4|≤
4114 m|β2|4α+↓|4|≤s|β2(|≤rX|β1|4α+↓|4{H11}|¬H{H10}|v41|4α_↓|4|≤r|
4114 g2|)X|β2),|E|;(25)|?{A9}|πthen |εY|β1 |πand |εY|β2
4120 |πare |εdependent |πrandom variables, normally
4125 distributed with menas |ε|≤m|β1, |≤m|β2 |πstandard
4131 deviations |ε|≤s|β1, |≤s|β2, |πand with correlation
4137 coe∃cient |ε|≤r. (|πFor a generalization to |εn
4144 |πvariables, see exercise 13.)|'{A12}|≡D|≡.|9|≡T|≡h|≡e
4149 |≡e|≡x|≡p|≡o|≡n|≡e|≡n|≡t|≡i|≡a|≡l |≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡
4150 t|≡i|≡o|≡n|≡.|9|4After uniform deviates and normal
4155 deviates, the next most important random quantity
4162 is an |εexponential deviate. |πSuch numbers occur
4169 in ``arrival time'' situations; for example,
4175 if a radioactive substance emits alpha particles
4182 at a rate so that one particle is emitted every
4192 |ε|≤m |πseconds on the average, then the time
4200 between two successive emissions has the exponential
4207 distribution with mean |ε|≤m. |π|πThis distribution
4213 is de_ned by the formula|'{A9}|εF(x)|4α=↓|41|4α_↓|4e|gα_↓|gx
4218 |g/|g|≤m,!!x|4|¬R|40.|E|;(26)|?{A9}{A3}|*/(|↔O)|9|\Logarithm
4221 method.|9|4|πClearly, if |εy|4α=↓|4F(x)|4α=↓|41|4α_↓|4e|gα_↓
4223 |gx|g/|g|≤m, x|4α=↓|4F|gα_↓|g1(y)|4α=↓|4|→α_↓|≤m|πln(|ε1|4α_
4224 ↓|4y). |πTherefore by Eq. (6), |→α_↓|ε|≤m|πln(|ε1|4α_↓|4U)
4230 |πhas the exponential distribution. Since |ε1|4α_↓|4U
4236 |πis uniformly distributed when |εU |πis, we
4243 conclude that|'{A9}|εX|4α=↓|4|→α_↓|≤m|πln|4|εU|E|;
4246 (27)|?{A9}|πis exponentially distributed with
4251 mean |ε|≤m. (|πThe case |εU|4α=↓|40 |πmust be
4258 avoided.)|'{A12}|ε|*/(|↔P)|9|\Random minimization
4261 method.|9|4|πWe saw in Algorithm F that there
4268 are simple and fast alternatives to calculating
4275 the logarithm of a uniform deviate. The following
4283 especially e∃cient approach has been developed
4289 by G. Marsaglia, M. Sibuya, and J. H. Ahrens.|'
4298 {A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡S (|εExponential
4301 distribution with mean |≤m).|9|4|πThis algorithm
4306 produces exponential deviates on a binary computer,
4313 using uniform deviates of |εm-|πbit accuracy.
4319 The constants|'{A9}|εQ[k]|4α=↓|4|(|πln|42|d21*3|)|4α+↓|4|((ln
4321 |42)|g2|d22*3|)|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4|((ln|42)|g|εk|d2k*3
4321 |)|4,!!k|4|¬R|41,|E|;(28)|?{A9}|πshould be precomputed,
4326 extending until |εQ[k]|4|¬Q|41|4α_↓|42|g1|gα_↓|gm.|'
4329 {A3}{I1.7H}|π|≡S|≡1|≡.|9[Get |εU |πand shift.]
4333 Generate a uniform random binary fraction |εU|4α=↓|4.b|β1b|β
4339 2|4.|4.|4.|4b|βm; |πLocate the _rst zero bit
4345 |εb|βj, |πand shift o= the leading |εj |πbits,
4353 setting |εU|5|¬L|5.b|βj|βα+↓|β1|4.|4.|4.|4b|βm.
4355 (|πAs in Algorithm F, the average value of |εj
4364 |πis 2.)|'{A3}|≡S|≡2|≡.|9[Immediate acceptance?]|9|4If
4368 |εU|4|¬W|4|πln|42, set |εX|5|¬L|5|ε|≤m(j|4|πln|42|4α+↓|4|εU)
4370 |πand terminate the algorithm. (Note that |εQ[1]|4α=↓|4|πln
4377 |42.)|'{A3}|≡S|≡3|≡.|9[Minimize.]|9|4Find the
4380 least |εk|4|¬R|42 |πsuch that |εU|4|¬W|4Q[k].
4385 |πGenerate |εk |πnew uniform deviates |εU|β1,|4.|4.|4.|4,|4U
4390 |βk |πand set |εV|5|¬L|5|πmin(|εU|β1,|4.|4.|4.|4,|4U|βk).|'
4394 {A3}|≡S|≡4|≡.|9[Deliver the answer.]|9|4|πSet
4397 |εX|5|¬L|5|≤v(j|4α+↓|4V)|πln|42.|'{A12}{IC}|≡E|≡.|9|≡O|≡t|≡h
4398 |≡e|≡r |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o|≡u|≡s |≡d|≡i|≡s|≡t|≡r|≡i|≡b|
4400 ≡u|≡t|≡i|≡o|≡n|≡s|≡.|9|4Let us now consider brie⊗y
4405 how to handle some other distributions which
4412 arise reasonably often in practice.|'{A12}|ε|\(|↔O)|9|\The
4418 gamma distribution |πof order |εa|4|¬Q|40 |πis
4424 de_ned by|'{A9}|εF(x)|4α=↓|4|(1|d2|≤)(a)|)|4|↔j|urx|)0|)|4t|
4426 ga|gα_↓|g1e|gα_↓|gtdt,!!x|4|¬R|40.|E|;(29)|?{A9}|πWhen
4429 |εa|4α=↓|41, |πthis is the exponential distribution
4435 with mean 1; when |εa|4α=↓|4|f1|d32|), |πit is
4442 the distribution of |f1|d32|)|εZ|g2, |πwhere
4447 |εZ |πhas the normal distribution (mena 0, variance
4455 1). If |εX |πand |εY |πare independent gamma-distributed
4463 random variables, of order |εa |πand |εb |πrespectively,
4471 then |εX|4α+↓|4Y |πhas the gamma distribution
4477 of order |εa|4α+↓|4b. |πThus, for example, the
4484 sum of |εk |πindependent exponential deviates
4490 with mean 1 has the gamma distribution of order
4499 |εk. |πIf the logarithm method (27) is being
4507 used to generate these exponential deviates,
4513 we need compute only one logarithm: |εX|4|¬L|4|→α_↓|πln(|εU|
4519 β1|4.|4.|4.|4U|βk), |πwhere |εU|β1,|4.|4.|4.|4,|4U|βk
4522 |πare nonzero uniform deviates. This technique
4528 handles all integer orders |εa; |πto complete
4535 the picture, a suitable method for 0|4|¬W|4a|4|¬W|41
4542 |πappears in exercise 16.|'!|9|7The simple logarithm
4549 method is much too slow when |εa |πis large,
4558 since it requires |"l|εa|"L |πuniform deviates.
4564 For large |εa, |πthe following remarkable algorithm
4571 due to J. H. Ahrens is reasonably e∃cient, and
4580 easy to write in terms of standard subroutines.|'
folio 165 galley 7
4588 |{U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addison-
4589 Wesley)|9f. 165|9ch. 3|9g. 7c|'{A6}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|
4593 ≡m |≡A (|εGamma distribution of order a|4|¬Q|41).|'
4600 {A3}{I1.8H}|π|≡A|≡1|≡.|9[Generate candidate.]
4602 Set |εY|4|¬L|4|πtan(|ε|≤pU), |πwhere |εU |πis
4607 a uniform deviate, and set |εX|5|¬L|5{H11}|¬H{H10}|v42a|4α_↓
4612 |41|)|2Y|4α+↓|4a|4α_↓|41. (|πIn place of tan(|ε|≤pU)
4617 |πwe could use a polar method, e.g., |εV|β2/V|β1
4625 |πin step P4 of Algorithm P.)|'{A3}|≡A|≡2|≡.|9[Accept?]
4632 If |εX|4|¬E|40, |πreturn to A1. Otherwise generate
4639 a new uniform deviate |εU, |πand return to A1
4648 if |εU|4|¬Q|4(1|4α+↓|4Y|g2)|πexp{H12}({H10}(|εa|4α_↓|41)|πln
4649 {H12}({H10}X/(a|4α_↓|41){H12}){H10}|4α_↓|4{H11}|¬H{H10}|v42a
4649 |4α_↓|41|)|2Y{H12}){H10}.|'{A12}{IC}|πThe average
4652 number of times step A1 is performed is |→|¬W1.902
4661 when |εa|4|¬R|43. |πFor further discussion, proof,
4667 and a slightly more complex method which is two
4676 to three times faster, see J. H. Ahrens and U.
4686 Dieter, |εComputing |π|≡1|≡2 (1974), 223<246.
4691 A still faster method, using normal and exponential
4699 deviates, has been developed by G. Marsaglia
4706 [to appear].|'{A12}|ε|*/(|"P)|9|\The beta distribution
4711 |πwith positive parameters |εa |πand |εb |πis
4718 de_ned by|'{A9}|εF(x)|4α=↓|4|(|≤)(a|4α+↓|4b)|d2|≤)(a)|≤)(b)|
4720 )|4|↔j|urx|)0|)|4t|ga|gα_↓|g1(1|4α_↓|4t)|gb|gα_↓|g1|4dt,!!0|
4720 4|¬E|4x|4|¬E|41.|E|;(30)|?{A9}|πLet |εX|β1 |πand
4725 |εX|β2 |πbe independent gamma deviates of order
4732 |εa |πand |εb, |πrespectively, and set |εX|5|¬L|5X|β1/(X|β1|
4738 4α+↓|4X|β2). |πAnother method, useful for small
4744 |εa |πand |εb, |πis to set |εY|β1|5|¬L|5U|ur1/a|)1|)
4751 |πand |εY|β2|5|¬L|5U|ur1/b|)2|) |πrepeatedly
4754 until |εY|β1|4α+↓|4Y|β2|4|¬E|41; |πthen |εX|5|¬L|5Y|β1/(Y|β1
4757 |4α+↓|4Y|β2). [|πSee M. D. J|=4ohnk, |εMetrika
4763 |π|≡8 (1964), 5<15.] |πStill another approach,
4769 if |εa |πand |εb |πare integers (not too large),
4778 is to set |εX |πto the |εb|πth largest of |εa|4α+↓|4b|4α_↓|4
4787 1 |πindependent uniform deviates (cf. exercise
4793 5<7).|'{A12}|ε|*/(|↔L)|9|\The chi-square distribution
4797 |πwith |ε|≤n |πdegrees of freedom (Eq. 3.3.1<22)
4804 is obtained by setting |εX|5|¬L|52Y, |πwhere
4810 |εY |πhas the gamma distribution of order |ε|≤n/2.|'
4818 {A12}|ε|*/(|↔C)|9|\The F-distribution (|πvariance-ratio
4821 distribution) with |εv|β1 |πand |εv|β2 |πdegrees
4827 of freedom, de_ned by|'{A9}|εF(x)|4α=↓|4|(|≤n|ur|≤n|β1/2|)1|
4831 )|≤n|ur|≤n|β2/2|)2|)|≤){H12}({H10}(|≤n|β1|4α+↓|4|≤n|β2)/2{H1
4831 2}){H10}|d2|≤)(|≤n|β1/2)|≤)(|≤n|β2/2)|)|4|↔j|urx|)0|)|4t|g|≤
4831 n|r1|g/|g2|gα_↓|g1(|≤n|β2|4α+↓|4|≤n|β1t)|gα_↓|g|≤n|r1|g/|g2|
4831 gα_↓|g|≤n|r2|g/|g2|4dt,|E|;(31)|?{A9}|πwhere
4834 |εx|4|¬R|40. |πLet |εY|β1 |πand |εY|β2 |πbe independent,
4841 having the chi-square distribution with |ε|≤n|β1,
4847 |≤n|β2 |πdegrees of freedom, respectively; set
4853 |εX|5|¬L|5Y|β1|≤n|β2/Y|β2|≤n|β1. |πOr set |εX|5|¬L|5|ε|≤n|β2
4856 Y/|≤n|β1(1|4α_↓|4Y), |πwhere |εY |πhas the beta
4862 distribution with parameters |ε|≤n|β1/2, |≤n|β2/2.|'
4867 {A12}|ε|*/(|↔C)|9|\The t-distribution |πwith |ε|≤n
4871 |πdegrees of freedom, de_ned by|'{A9}|εF(x)|4α=↓|4|(|≤){H12}
4876 ({H10}(|≤n|4α+↓|41)/2{H12}){H10}|d2{H11}|¬H{H10}|v2|≤p|≤n|)|
4876 4|≤)(|≤n/2)|)|4|↔j|urx|)α_↓|¬X|)|4(1|4α+↓|4t|g2/|≤n)|gα_↓|g(
4876 |g|≤n|gα+↓|g1|g)|g/|g2|4dt.|E|;(32)|?{A9}|πLet
4879 |εY|β1 |πbe a nromal deviate (mean 0, variance
4887 1) and let |εY|β2 |πbe independent of |εY|β1,
4895 |πhaving the chi-square distribution with |ε|≤n
4901 |πdegrees of freedom; set |εX|5|¬L|5Y|β1/{H11}|¬H{H10}|v4Y|β
4905 2/|≤n|).|'{A12}|*/(|↔o)|9|\Random point on n-dimensional
4910 sphere with radius one.|9|4|πLet X|β1, X|β2,|4.|4.|4.|4,
4916 X|βn |πbe independent normal deviates (mean 0,
4923 variance 1); the desired point on the unit sphere
4932 is|'{A9}|ε(X|β1/r,|6X|β2/r,|4.|4.|4.|4,|6X|βn/r),!!|πwhere!!
4933 |εr|4α=↓|4{H11}|¬H{H10}|v4X|=|β1|g2|4α+↓|4X|=|β2|g2|4α+↓|4|¬
4933 O|4|¬O|4|¬O|4α+↓|4X|=|βn|g2|).|E|;(33)|?{A9}|πNote
4936 that if the |εX'|πs are calculated using the
4944 polar method, Algorithm P, we compute two independent
4952 |εX'|πs each time, and |εX|=|β1|g2|4α+↓|4X|=|β2|g2|4α=↓|4|→α
4956 _↓2|4|πln|4|εS (|πin the notation of that algorithm);
4963 this saves a little of the time needed to evaluate
4973 |εr. |πThe validity of this method comes from
4981 the fact that the distribution function for the
4989 point (|εX|β1,|4.|4.|4.|4, X|βn) |πhas a density
4995 which depends only on the distance from the origin,
5004 so when it is projected onto the unit sphere
5013 it has the uniform distribution. this method
5020 was _rst suggested by G. W. Brown, in |εModern
5029 Mathematics for the Engineer, |πFirst series,
5035 ed. by E. F. Beckenbach (New York: McGraw-Hill,
5043 1956), p. 302. To get a random point |εinside
5052 |πthe |εn-|πsphere, R. P. Brent suggests taking
5059 a point on the surface and multiplying it by
5068 |εU|g1|g/|gn.|'|π!|9|7In three dimensions a signi_cantly
5074 simpler method can be used, since each individual
5082 coordinate is uniformly distributed between |→α_↓1
5088 and 1: Find |εV|β1, V|β2, S |πby steps P1<P3
5097 of Algorithm P; then the desired point is (|ε|≤aV|β1,
5106 |≤aV|β2, 2S-1), |πwhere |ε|≤a|4α=↓|42{H11}|¬H{H10}|v41|4α_↓|
5109 4S|). [|πRobert E. Knop, |εCACM |π|≡1|≡3 (1970),
5116 326.]|'{A12}|≡F|≡.|9|≡I|≡m|≡p|≡o|≡r|≡t|≡a|≡n|≡t
5118 |≡i|≡n|≡t|≡e|≡g|≡e|≡r|≡-|≡v|≡a|≡l|≡u|≡e|≡d |≡d|≡i|≡s|≡t|≡r|≡
5119 i|≡b|≡u|≡t|≡i|≡o|≡n|≡s|≡.|9|4A probability distribution
5122 which consists only of integer values can essentially
5130 be handled by the techniques described at the
5138 beginning of this section; but some of these
5146 distributions are so important in practice, they
5153 deserve special mention here.|'{A12}|\|ε(|↔O)|9|\The
5158 geometric distribution.|9|4|πIf some event occurs
5163 with probability |εp, |πthe number |εN |πof independent
5171 trials needed until the event _rst occurs (or
5179 between occurrences of the event) has the geometric
5187 distribution. We have |εN|4α=↓|41 |πwith probability
5193 |εp, N|4α=↓|42 |πwith probability (|ε1|4α_↓|4p)p,|4.|4.|4.|4
5197 , N|4α=↓|4n |πwith probability (|ε1|4α_↓|4p)|gn|gα_↓|g1p.
5202 (|πThis is essentially the same situation as
5209 we have already considered in the gap test of
5218 Section 3.3.2; it is also directly related to
5226 the number of times certain loops in the algorithms
5235 of this section are executed, e.g., steps P1<P3
5243 of the polar method.)|'!|9|7A convenient way
5250 to generate a variable with this distribution
5257 is to set|'{A9}|εN|5|¬L|5|"p|πln|4|εU/|πln(|ε1|4α_↓|4p)|"P.|
5260 E|;(34)|?{A9}|πTo check this formula, we observe
5268 that |"pln|4|εU/|πln(|ε1|4α_↓|4p)|"P|4α=↓|4n
5270 |πif and only if |εn|4α_↓|41|4|¬W|4|πln|4|εU/|πln(1|4α_↓|4|ε
5274 p)|4|¬E|4n, |πthat is, (|ε1|4α_↓|4p)|gn|gα_↓|g1|4|¬Q|4U|4|¬R
5277 |4(1|4α_↓|4p)|gn, |πand this happens with probability
5283 |εp(1|4α_↓|4p)|gn|gα_↓|g1 |πas required. Note
5287 thatln|4|εU |πcan be replaced by |→α_↓|εY, |πwhere
5294 |εY |πhas the exponential distribution with mean
5301 1.|'!|9|7The special case |εp|4α=↓|4|f1|d32|)
5306 can be handled more easily on a binary computer,
5315 since formula (34) becomes |εN|4|¬L|4|"p|→α_↓|πlog|β2|4|εU|¬
5319 P, |πthat is, |εN |πis one more than the number
5329 of leading zero bits in the binary representation
5337 of {U0}{H10L12M29}W58320#Computer Programming|9(Knuth/Addiso
folio 173 galley 8
5339 n-Wesley)|9f. 173|9ch. 3|9g. 8c|'{A6}|ε|*/(|↔P)|9|\The
5344 binomial distribution (t,|4p).|9|4|πIf some event
5349 occurs with probability |εP |πand if we carry
5357 out |εt |πindependent trials, the total number
5364 |εN |πof occurrences equals |εn |πwith probability
5371 (|ε|=|βn|gt)p|gn(1|4α_↓|4p)|gt|gα_↓|gn. (|πSee
5373 Section 1.2.10.) |πIn other words if we generate
5381 |εU|β1,|4.|4.|4.|4, U|βt, |πwe want to count
5387 how many of these are |→|¬W|εp. |πFor small |εt
5396 |πwe can obtain |εN |πin exactly this way.|'!|9|7For
5405 large |εt, |πwe can generate a beta variate |εX
5414 |πwith integer parameters |εa |πand |εb |πwhere
5421 |εa|4α+↓|4b|4α_↓|41|4α=↓|4t; |πthis e=ectively
5424 gives us the |εb|πth largest of |εt |πelements,
5432 without bothering to generate the other elements.
5439 Now if |εX|4|¬R|4p, |πwe set |εN|4|¬L|4N|β1 |πwhere
5446 |εN|β1 |πhas the binomial distribution (|εa|4α_↓|41,|4p/X),
5452 |πsince this tells us how many of |εa|4α_↓|41
5460 |πrandom numbers in the range [0,|4|εX) |πare
5467 |→|¬W|εp; |πand if |εX|4|¬W|4p, |πwe set |εN|4|¬L|4a|4α+↓|4N
5473 |β1 |πwhere |εN|β1 |πhas the binomial distribution
5480 {H12}({H10}|εb|4α_↓|41,|4(p|4α_↓|4X)/(1|4α_↓|4X){H12}){H10},
5480 |πsince |εN|β1 |πtells us how many of |εb|4α_↓|41
5489 |πrandom numbers in the range [|εX,|41) |πare
5496 |ε|→|¬Wp. |πBy choosing |εa|4α=↓|41|4α+↓|4|"lT/2|"L,
5500 |πthe parameter |εt |πwill be reduced to a reasonable
5509 size after about lg|4|εt |πreductions of this
5516 kind. (This approach is due to J. H. Ahrens,
5525 who has also suggested an alternative procedure
5532 for medium-sized |εt; |πsee exercise 27.)|'{A12}|ε|*/(|↔L)|9|
5538 \The Poisson distribution |πwith mean |ε|≤m.|9|4|πThis
5544 distribution is related to the exponential distribution
5551 as the binomial distribution is related to the
5559 geometric: it represents the number of occurrences,
5566 per unit time, of an event which can occur at
5576 any instant of time; for example, the number
5584 of alpha particles emitted by radioactive substance
5591 in a single second has a Poisson distribution.|'
5599 !|9|7|πAccording to this principle, we can produce
5606 a Poisson deviate |εN |πby generating independent
5613 exponential deviates |εX|β1,|4X|β2,|4.|4.|4.
5616 |πwith mean 1/|ε|≤m, |πstopping as soon as |εX|β1|4α+↓|4|¬O|
5623 4|¬O|4|¬O|4α+↓|4X|βm|4|¬R|41; |πthen |εN|4|¬L|4m|4α_↓|41.
5626 |πThe probability that |εX|β1|4α+↓|4|¬O|4|¬O|4|¬O|4X|βm|4|¬R
5629 |41 |πis the probability that a gamma deviate
5637 of order |εm |πis |ε|→|¬R|≤m, |πnamely |ε|↔j|ur|¬X|)|≤m|)|4t
5643 |gm|gα_↓|g1e|gα_↓|gt|4dt/(m|4α_↓|41)*3; |πhence
5645 the probability that |εN|4α=↓|4n |πis|'{A9}|ε|(1|d2n*3|)|4|↔j
5650 |ur|¬X|)|≤m|)|4t|gne|gα_↓|gt|4dt|4α_↓|4|(1|d2(n|4α_↓|41)*3|)|
5650 4|↔j|ur|¬X|)|≤m|)|4t|gn|gα_↓|g1e|gα_↓|gt|4dt|4α=↓|4e|gα_↓|g|
5650 ≤m|4|(|≤m|gn|d2n*3|)|4,!!n|4|¬R|40.|E|'(35)|?{A9}|πIf
5653 we generate exponential deviates by the logarithm
5660 method, the above recipe tells us to stop when
5669 |→α_↓(ln|4|εU|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4|πln|4|εU|βm)/|≤m
5669 |4|¬R|41. |πSimplifying this expression, we see
5675 that the desired Poisson deviate can be obtained
5683 by calculating |εe|gα_↓|g|≤m, |πconverting it
5688 to a _xed-point representation, then generating
5694 one or more uniform deviates |εU|β1,|4U|β2,|4.|4.|4.
5700 |πuntil the product |εU|β1|4.|4.|4.|4U|βm|4|¬E|4e|gα_↓|g|≤m,
5703 |πand _nally setting |εN|4|¬L|4m|4α_↓|41. |πOn
5709 the average this requires the generation of |ε|≤m|4α+↓|41
5717 |πuniform deviates, so it is a very useful approach
5726 when |ε|≤m |πis not too large. For very small
5735 |ε|≤m |πit would be better to generate exponential
5743 variables directly with a streamlined version
5749 of Algorithm S.|'!|9|7When |ε|≤m |πis large,
5756 we can obtain a method of order log|4|ε|≤m |πby
5765 using the fact that we know how to handle the
5775 gamma and binomial distributions for large orders:
5782 First generate |εX |πwith the gamma distribution
5789 of order |εm|4α=↓|4|"l|≤a|≤m|"L, |πwhere |ε|≤a
5794 |πis a suitable constant. (Since |εX |πis equivalent
5802 to |→α_↓ln(|εU|β1|4.|4.|4.|4U|βm), |πwe are essentially
5807 bypassing |εm |πsteps of the previous method.)
5814 If |εX|4|¬W|4|≤m, |πset |εN|4|¬L|4m|4α+↓|4N|β1,
5818 |πwhere |εN|β1 |πis a Poisson deviate of mean
5826 |ε|≤m|4α_↓|4X; |πand if |εX|4|¬R|4|≤m, |πset
5831 |εN|4|¬L|4N|β1, |πwhere |εN|β1 |πhas the binomial
5837 distribution (|εm|4α_↓|41,|4|≤m/X). |πThis method
5841 is due to J. H. Ahrens and U. Dieter, whose experiments
5852 suggest that |f7|d38|) is a good choice for |ε|≤a.|'
5861 !|9|7|πThe validity of the above reduction when
5868 |εX|4|¬R|4|≤m |πis a consequence of the following
5875 important principle: ``Let |εX|β1,|4.|4.|4.|4,
5879 X|βm |πbe independent exponential deviates with
5885 the same mean; let |εS|βj|4α=↓|4X|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+
5889 ↓|4X|βj |πand let |εV|βj|4α=↓|4S|βj/S|βm |πfor
5894 |ε1|4|¬E|4j|4|¬E|4m. |πThen the distribution
5898 of |εV|β1,|4V|β2,|4.|4.|4.|4, V|βm|βα_↓|β1 |πis
5902 the same as the distribution of |εm|4α_↓|41 |πindependent
5910 uniform deviates sorted into increasing order.''
5916 To establish this principle formally, we compute
5923 the probability that |εV|β1|4|¬E|4v|β1,|4.|4.|4.|4,
5927 V|βm|βα_↓|β1|4|¬E|4v|βm|βα_↓|β1, |πgiven the
5930 value of |εS|βm|4α=↓|4s, |πfor arbitrarily values
5936 |ε0|4|¬E|4v|β1|4|¬E|4|¬O|4|¬O|4|¬O|4|¬Ev|βm|βα_↓|β1|4|¬E|41:
5936 |'{A9}{M36}|↔j|urv|β1s|)0|)|4|≤m|πe|ε|gα_↓|gt|r1|g/|g|≤m|4dt
5937 |β1|4|↔j|urv|β2sα_↓t|β1|)0|)|3|≤me|gα_↓|gt|r2|g/|g|≤m|4dt|β2
5937 |4.|4.|4.|4|↔j|urv|βm|βα_↓|β1sα_↓t|β1|4α_↓|4|¬O|4|¬O|4|¬O|4α
5937 _↓|4t|βm|βα_↓|β2|3|≤me|gα_↓|gt|rm|rα_↓|r1|g/|g|≤m|4dt|βm|βα_
5937 ↓|β1|4|¬O|4|≤me|gα_↓|g(|gs|gα_↓|gt|r1|gα_↓|4|g|¬O|4|g|¬O|4|g
5937 |¬O|4|gα_↓|gt|rm|rα_↓|r1|g)|g/|g|≤m|d2|↔j|urs|)0|)|3|ε|≤me|g
5937 α_↓|gt|r1|g/|g|≤m|4dt|β1|4|↔j|ursα_↓t|β1|)0|)|3|≤me|gα_↓|gt|
5937 rz|g/|g|≤m|4dt|β2|4.|4.|4.|3|↔j|ursα_↓t|β1α_↓|4|¬O|4|¬O|4|¬O
5937 |4α_↓t|βm|βα_↓|β2|)0|)|3|≤me|gα_↓|gt|rm|rα_↓|r1|g/|g|≤m|4dt|
5937 βm|βα_↓|β1|4|¬O|4|≤me|gα_↓|g(|gs|gα_↓|gt|r1|gα_↓|g|4|g|¬O|g|
5937 4|g|¬O|g|4|g|¬O|g|4|gα_↓|gt|rm|rα_↓|r1|g)|g/|g|≤m|)|'
5938 {A6}{M29}α=↓|4|(|↔j|urv|β1|)0|)|3du|β1|3|↔j|urv|β2|)u|β1|)|3
5938 du|β2|4.|4.|4.|3|↔j|urv|βm|βα_↓|β1|)u|βm|βα_↓|β2|)|3du|βm|βα
5938 _↓|β1|d2|↔j|ur1|)0|)|3|↔j|ur1|)u|β1|)|3du|β2|4.|4.|4.|3|↔j|u
5938 r1|)u|βm|βα_↓|β2|)|3du|βm|βα_↓|β1|)|4,|?{A9}|πby
5940 making the substitution |εt|β1|4α=↓|4su|β1, t|β1|4α+↓|4t|β2|
5944 4α=↓|4su|β2,|4.|4.|4.|4,|4t|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4t|β
5944 m|βα_↓|β1|4α=↓|4su|βm|βα_↓|β1. |πThe latter ratio
5948 is the corresponding probability that uniform
5954 deviates |εU|β1,|4.|4.|4.|4, U|βm|βα_↓|β1 |πsatisfy
5958 |εU|β1|4|¬E|4v|β1,|4.|4.|4.|4, U|βm|βα_↓|β1|4|¬E|4v|βm|βα_↓|
5959 β1, |πgiven that |εU|β1|4|¬E|4|¬O|4|¬O|4|¬O|4|¬E|4U|βm|βα_↓|
5962 β1.|'!|9|7|πA more e∃cient but somewhat more
5969 complex technique for binomial and Poisson deviates
5976 is sketched in exercise 22.|'{A12}|≡G|≡.|9|≡F|≡o|≡r
5982 |≡f|≡u|≡r|≡t|≡h|≡e|≡r |≡r|≡e|≡a|≡d|≡i|≡n|≡g|≡.|9|4The
5984 book |εNon-Uniform Random Numbers |πby J. H.
5991 Ahrens and U. Dieter (New York: Wiley 1976) discusses
6000 many more algorithms for the generation of random
6008 variables with nonuniform distributions, together
6013 with a careful consideration of the e∃ciency
6020 of each technique on typical computers.|'!|9|7Exercise
6027 16 is recommended as a review of many of the
6037 techniques in this section.|'|Ha*?*?{U0}{H10L12M29}W58320#Comp
folio 176 galley 9
6041 uter Programming|9(Knuth/Addison-Wesley)|9f.
6043 176|9ch. 3|9g. 9c|'{A6}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'
6047 {A12}{H9L11}|5|5|≡1|≡.|9|4[|ε|*/|↔O|↔c|\] |πIf
6049 |ε|≤a |πand |ε|≤b |πare real numbers with |ε|≤a|4|¬W|4|≤b,
6057 |πhow would you generate a random real number
6065 uniformly distributed between |ε|≤a |πand |ε|≤b?|'
6071 {A3}|5|5|≡2|≡.|9|4[|*/M16|\] |πAssuming that |εmU
6075 |πis a random integer between 0 and |εm|4α_↓|41,
6083 |πwhat is the |εexact |πprobability that |ε|"lkU|"L|4α=↓|4r,
6089 |πif 0|4|¬E|4|εr|4|¬W|4k? |πCompare this with
6095 the desired probability |ε1/k.|'{A3}|5|5|≡3|≡.|9|4|ε[|*/|↔O|↔
6099 M|\] |πDiscuss treating |εU |πas an integer and
6107 |εdividing |πby |εk |πto get a random integer
6115 between 0 and |εk|4α_↓|41, |πinstead of multiplying
6122 as suggested in the text. Thus (1) would be changed
6132 to|'{A9}|¬e|¬n|¬t|¬a!|¬o|;|¬l|¬d|¬x|6|6|¬u|;|¬d|¬i|¬v|6|6|¬k
6135 |;{A6}with the result appearing in register X.
6143 Is this a good method?|'{A3}|5|5|≡4|≡.|9|4[|ε|*/M|↔P|↔c|\]
6149 |πProve the two relations in (7).|'{A3}|5|5|≡5|≡.|9|4[|ε|*/21
6155 |\] |πSuggest an e∃cient method to compute a
6163 random variable with the distribution |εpx|4α+↓|4qx|g2|4α+↓|
6168 4rx|g3, |πwhere |εp|4|¬R|40, q|4|¬R|40, r|4|¬R|40,
6173 |πand |εp|4α+↓|4q|4α+↓|4r|4α=↓|41.|'{A3}|5|56.|9|4[|ε|*/HM|↔P
6175 |↔O|\] |πA quantity |εX |πis computed by the
6183 following method:|'{A9}!!|2``|εStep 1.|9|4|πGenerate
6187 two independent uniform deviates |εU, V.|'{A2}!!|9|3Step
6194 2.|9|4|πIf |εU|g2|4α+↓|4V|g2|4|¬R|41, |πreturn
6197 to step 1; otherwise set |εX|4|¬L|4U.''|'{A9}|πWhat
6204 is the distribution function of |εX? |πHow many
6212 times will step 1 be performed? (Give the mean
6221 and standard deviation.)|'{A3}|5|5|≡7|≡.|9|4[|ε|*/M|↔O|↔l|\]
6225 |πExplain why, in Marsaglia's method for normal
6232 deviates, the decision to make |εp|βj |πa multiple
6240 of |f1|d3256|) implies that we should take |εp|βj|4α=↓|4|"l6
6247 4f|1(j/4)|"L/256, 1|4|¬E|4j|4|¬E|412.|'{A3}|5|5|≡8|≡.|9|4[|*/
6249 |↔O|↔c|\] |πWhy are there ``skinny'' rectangles
6255 |εf|β1|β3,|4.|4.|4.|4,|4f|β2|β4 |πin Marsaglia's
6258 method, as well as the larger rectangles |εf|β1,|4.|4.|4.|4,
6265 |4f|β1|β2? (|πWhy is this better than putting
6272 (|εf|β1,|4f|β1|β3), (f|β2,|4f|β1|β4),|4.|4.|4.
6274 |πtogether into single large rectangles?)|'{A3}|5|5|≡9|≡.|9|
6279 4[|ε|*/HM|↔O|↔c|\] |πWhy is the curve |εf|1(x)
6285 |πof Fig. 9 concave downward for |εx|4|¬W|41,
6292 |πconcave upward for |εx|4|¬Q|41?|'{A3}|≡1|≡0|≡.|9|4[|ε|*/HM|
6296 ↔P|↔L|\] |πGive formulas for the numbers |εa|βj,
6303 b|βj, p|βj|βα+↓|β2|β4 |πsuch that the densities
6309 |εf|βj|βα+↓|β2|β4(x) |πof Fig. 9 can be calculated
6316 by Algorithm L with |εa|4α=↓|4a|βj, b|4α=↓|4b|βj,
6322 h|4α=↓|4|f1|d34|), |πand |εs|4α=↓|4|f1|d34|)(j|4α_↓|41),
6325 |πfor |ε1|4|¬E|4j|4|¬E|412.|'{A3}|≡1|≡1|≡.|9|4[|*/HM|↔P|↔p|\]
6327 |πProve that steps M8<M9 of Algorithm M generate
6336 a random variable with the appropriate tail of
6344 the normal distribution; i.e., the probability
6350 that |εX|4|¬E|4x |πshould be|'{A9}|ε|↔j|urx|)3|)|3e|gα_↓|gt|
6354 i2|g/|g2|4dt|↔h|2|↔j|ur|¬X|)3|)|3e|gα_↓|gt|i2|g/|g2|4dt,!!x|
6354 4|¬R|43.|;{A9}|πwhen |εx|4|¬R|43.|'{A3}|≡1|≡2|≡.|9|4[|*/HM|↔M
6357 |↔o|\] |πTeichroew's polynomial, Eq. (22), is
6363 a truncated sum of Chebyshev polynomials and
6370 it may not be the best possible polynomial approximation
6379 of degree 9|'{A12}[|εHint|*/:|9|\|πIt is a special
6386 case of the rejection method, with |εg(t)|4α=↓|4Cte|gα_↓|gt|
6392 i2|g/|g2 |πfor some |εC.]|'{A3}|≡1|≡2|≡.|9|4[|ε|*/HM|↔P|↔L|\]
6396 |π(R. P. Brent.) Prove that the numbers |εa|βj
6405 |πde_ned in (23) satisfy |εa|=|βj|g2/2|4α_↓|4a|ur2|)jα_↓1|)/
6409 2|4|¬W|4|πln|42 for all |εj|4|¬R|41. [Hint|*/:
6414 |\|πIf |εf|1(x)|4α=↓|4e|gx|i2|g/|g2|4|¬J|ur|¬X|)x|)|4e|gα_↓|
6415 gt|i2|g/|g2|4dt, |πshow that |εf|1(x)|4|¬W|4f|1(y)
6419 |πfor |ε0|4|¬E|4x|4|¬W|4y.]|'{A3}|≡1|≡3|≡.|9|4[|ε|*/HM|↔P|↔C|
6421 \] |πGiven a set of |εn |πindependent normal
6429 deviates, |εX|β1, X|β2,|4.|4.|4.|4,|4X|βn, |πwith
6433 mean 0 and variance 1, show how to _nd constants
6443 |εb|βj |πand |εa|βi|βj, 1|4|¬E|4j|4|¬E|4i|4|¬E|4n,
6447 |πso that if|'{A9}!!|εY|β1|4α=↓|4b|β1|4α+↓|4a|β1|β1X|β1,!!Y|
6450 β2|4α=↓|4b|β2|4α+↓|4a|β2|β1X|β1|4α+↓|4a|β2|β2X|β2,!!.|4.|4.|
6450 4,|'{A4}|εY|βn|4α=↓|4b|βn|4α+↓|4a|βn|β1X|β1|4α+↓|4a|βn|β2X|β
6451 2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓a|βn|βnX|βn,!!|?{A9}|πthen
6453 |εY|β1, Y|β2,|4.|4.|4.|4,|4Y|βn |πare dependent,
6457 normally distributed variables, |εY|βj |πhas
6462 mean |ε|≤m|βj, |πand the |εY'|πs have a given
6470 covariance matrix (|εc|βi|βj). (|πThe covariance,
6475 |εc|βi|βj, |πof |εY|βi |πand |εY|βj |πis de_ned
6482 to be the average value of (|εY|βi|4α_↓|4|≤m|βi)(Y|βj|4α_↓|4
6488 |≤m|βj). |πIn particular, |εc|βj|βj |πis the
6494 variance of |εY|βj, |πthe square of its standard
6502 deviation. Not all matrices (|εc|βi|βj) |πcan
6508 be covariance matrices, and your construction
6514 is, of course, only supposed to work whenever
6522 a solution to the given conditions is possible.)|'
6530 {A3}|≡1|≡4|≡.|9|4[|ε|*/M|↔P|↔O|\] |πIf |εX |πis
6534 a random variable with continuous distribution
6540 |εF(x), |πand if |εc |πis a constant, what is
6549 the distribution of |εcX?|'{A3}|≡1|≡5|≡.|9|4[|*/HM21|\]
6554 |πIf |εX|β1 |πand |εX|β2 |πare independent random
6561 variables with the respective distributions |εF|β1(x)
6567 |πand |εF|β2(x), |πand the densities |εf|β1(x)|4α=↓|4F|=|β1|
6572 ¬S(x), f|β2(x)|4α=↓|4F|=|β2|¬S(x), |πwhat are
6576 the distribution and density functions of the
6583 quantity |εX|β1|4α+↓|4X|β2?|'{A3}|≡1|≡6|≡.|9|4[|ε|*/HM22|\]
6586 |π(J. H. Ahrens.) Develop an algorithm for gamma
6594 deviates of order |εa |πwhen |ε0|4|¬W|4|εa|4|¬E|41,
6600 |πusing the rejection method with |εcg(t)|4α=↓|4t|ga|gα_↓|g1
6605 /|≤)(a) |πfor |ε0|4|¬W|4t|4|¬W|41, |εe|gα_↓|gt/|≤)(a)
6609 |πfor |εt|4|¬R|41.|'{A3}|≡1|≡7|≡.|9|4[|ε|*/M|↔P|↔M|\]
6612 |πWhat is the |εdistribution function F(x) |πfor
6619 the geometric distribution with probability |εp?
6625 |πWhat is the |εgenerating function G(z)? |πWhat
6632 are the mean and standard deviation of this distribution?|'
6641 {A3}|≡1|≡8|≡.!|9|7[|ε|*/N|↔P|↔M|\] |πSuggest a
6644 method to compute a random integer |εN |πfor
6652 which |εN |πtakes the value |εn |πwith probability
6660 |εnp|g2(1|4α_↓|4p)|gn|gα_↓|g1, n|4|¬R|40. (|πThe
6663 case of particular interest is when |εp |πis
6671 rather small.)|'{A3}|≡1|≡9|≡.|9|4[|ε|*/|↔P|↔P|\]
6674 |πThe |εnegative binomial distribution (t,|4p)
6679 |πhas integer values |εN|4α=↓|4n |πwith probability
6685 (|ε|ft|2α_↓|21|2α+↓|2n|d5n|))p|gt(1|4α_↓|4p)|gn.
6686 (|πUnlike the ordinary binomial distribution,
6691 |εt |πneed not be an integer, since this quantity
6700 is nonnegative for all |εn |πwhenever |εt|4|¬Q|40.)
6707 |πGeneralizing exercise 18, explain how to generate
6714 integers |εN |πwith this distribution when |εt
6721 |πis a small positive integer. What method would
6729 you suggest if |εt|4α=↓|4p|4α=↓|4|f1|d32|)?|'
6733 {A3}|≡2|≡0|≡.|9|4[|ε|*/|↔P|↔L|\] |π(``Marsaglia
6735 tables.'') Suppose we want to compute a random
6743 variable |εX |πwith the following _nite distribution:|'
6750 {A9}|hWith|6probability|4|∂α=↓|4|f1|d3512|)!|f8|d3512|)!|f3|
6750 d3512|)!|f8|d3512|)!|f3|d3512|)!|f6|d3512|)|E|n|;
6751 | |πValue|6of|6|εX|4|Lα=↓|4|5x|β1|6!|5x|β2|6!|5x|β3|6!|5x|β4
6751 |6!|5x|β5|6!|5x|β6>{A4}| |πWith|6probability|4|Lα=↓|4|f90|d3
6752 512|)!|f81|d3512|)!|f131|d3512|)!|f10|d3512|)!|f32|d3512|)!|
6752 f168|d3512|)>{A9}|πShow how the value of |εX
6759 |πcan be e∃ciently obtained from a random nine-bit
6767 binary integer, |εb|β1b|β2b|β3b|β4b|β5b|β6b|β7b|β8b|β9,
6770 |πusing a method analogous to (16) and using
6778 auxiliary tables containing less than 35 entries
6785 in all.|'|Ha*?*?*?*?{U0}{H10L12M29}W58320#Computer
folio 178 galley 10
6788 Programming|9(Knuth/Addison-Wesley)|9f. 178|9ch.
6790 3|9g. 10c|'{A6}{H9L11}|π|≡2|≡1|≡.|9|4[|ε|*/|↔P|↔M|\]
6793 |πThe situation in the preceding exercise is
6800 idealized in the sense that we usually will not
6809 have all probabilities an exact multiple of a
6817 convenient number like |f1|d3512|). Suppose the
6823 probabilities which appear in that exercise are
6830 modi_ed so that they are actually equal to the
6839 following:|'{A9}{A16}{M32}With|4probability|4|∂α=↓|4|f89|d35
6840 12|4α+↓|4|≤e|β1!|∂|f80|d3512|)|4α+↓|4|≤e|β2!|∂|f131|d3512|)|
6840 4α+↓|4|≤e|β3!|∂|f10|d3512|)|4α+↓|4|≤e|β4!|∂|f32|d3512|)|4α+↓
6840 |4|≤e|β5!|∂|f168|d3512|)|4α+↓|4|≤e|β6|E|'{B16}| Value|4of|4|
6841 εX|4|Lα=↓!|9x|β1|L!|9x|β2|L!|9x|β3|L!|9x|β4|L!|9x|β5|L!|9x|β
6841 6>{A25}{M29}|πHere |ε0|4|¬E|4|≤e|βj|4|¬W|4|f1|d3512|),
6844 |πand |ε|≤e|β1|4α+↓|4|≤e|β2|4α+↓|4|≤e|β3|4α+↓|4|≤e|β4|4α+↓|4
6845 |≤e|β5|4α+↓|4|≤e|β6|4α=↓|4|f2|d3512|). |πShow
6847 how to make this choice e∃ciently, when given
6855 a uniform binary deviate |εU|4α=↓|4.|εb|β1b|β2|4.|4.|4.|4b|β
6859 m |πwhere |εm |πis relatively large. As in exercise
6868 20, try to do this using auxiliary tables containing
6877 less than 35 entries in all.|'{A3}|≡2|≡2|≡.|9|4[|ε|*/HM|↔P|↔P
6883 |\] |πCan the exact Poisson distribution for
6890 large |ε|≤m |πbe obtained by generating an appropriate
6898 normal deviate, converting it to an integer in
6906 some convenient way, and applying a (possibly
6913 complicated) correction a small percent of the
6920 time?|'{A3}|≡2|≡3|≡.|9|4[|ε|*/HM|↔P|↔L|\] (|πJ.
6923 von Neumann.) Are the following two ways to generate
6932 a random quantity |εX |πequivalent (i.e., does
6939 the quantity |εX |πhave the same distribution)?|'
6946 {A3}!!|2|εMethod|6|ε|*/|↔O|\:|5|5|πSet |εX|5|¬L|5|πsin{H12}({
6947 H10}(|ε|≤p/2)U{H12}){H10}, |πwhere |εU |πis uniform.|'
6952 {A3}{I1.7l}|εMethod|6|*/|↔P|\:|5|5|πGenerate two
6954 uniform deviates, |εU, V, |πand if |εU|g2|4α+↓|4V|g2|4|¬R|41
6960 , |πrepeat until |εU|g2|4α+↓|4V|g2|4|¬W|41. |πThen
6965 set |εX|5|¬L|5|¬GU|g2|4α_↓|4V|g2|¬G/(U|g2|4α+↓|4V|β2).|'
6967 {A3}{IC}|π|≡2|≡4|≡.|9|4[|ε|*/HM|↔M|↔c|\] (|πS.
6969 Ulam, J. von Neumann.) Let |εV|β0 |πbe a randomly
6978 selected real number between 0 and 1, and de_ne
6987 the sequence |¬D|εV|βn|¬F |πby the rule |εV|βn|βα+↓|β1|4α=↓|
6993 44V|βn(1|4α_↓|4V|βn). |πNow if this computation
6998 is done with perfect accuracy, the result should
7006 be a sequence with the distribution sin|g2|4|ε|≤pU,
7013 |πwhere |εU |πis uniform, i.e., with distribution
7020 function|'{A9}|εF(x)|4α=↓|3|↔j|urx|)0|)|3dx/{H11}|¬H{H10}|v2
7021 2|≤px(1|4α_↓|4x)|).|;{A9}|πFor if we write |εV|βn|4α=↓|4|πsi
7026 n|g2|4|ε|≤pU|βn, |πwe _nd that |εU|βn|βα+↓|β1|4α=↓|4(2U|βn)|
7030 πmod 1; and by the fact that almost all real
7040 numbers have a random binary expansion (see Section
7048 3.5), this sequence |εU|βn |πis equidistributed.
7054 But if the computation of |εV|βn |πis done with
7063 only _nite accuracy, the above argument breaks
7070 down because we soon are dealing with noise from
7079 the roundo= error. [|εReference|*/: |\|πvon Neumann's
7085 |εCollected Works |π|≡5|≡, 768<770.]|'!|9|7Analyze
7090 the sequence |¬D|εV|βn|¬F |πde_ned above when
7096 only _nite accuracy is present, both empirically
7103 (for various di=erent choices of |εV|β0) |πand
7110 theoretically. Does the sequence have a distribution
7117 resembling the expected distribution?|'{A3}|≡2|≡5|≡.|9|4[|ε|
7121 */M|↔P|↔C|\] |πLet |εX|β1, X|β2,|4.|4.|4.|4,|4X|β5
7125 |πbe binary words each of whose bits is independently
7134 0 or 1 with probability |f1|d32|). What is the
7143 probability that a given bit position of |εX|β1|4|↔I|4{H12}(
7150 {H10}X|β2|4|↔i|4{H12}({H10}X|β3|4|↔I|4(X|β4|4|↔i|4X|β5){H12}
7150 )){H10} |πcontains a 1? Generalize.|'{U31}1|'
7156 {H9L11M29}{A3}|π|≡2|≡6|≡.|9|4[|ε|*/M|↔O|↔l|\]
7157 |πLet |εN|β1 |πand |εN|β2 |πbe independent Poisson
7164 deviates with respective means |ε|≤m|β1 |πand
7170 |ε|≤m|β2, |πwhere |ε|≤m|β1|4|¬Q|4|≤m|β2|4|¬R|40.
7173 |πProve or disprove: (a) |εN|β1|4α+↓|4N|β2 |πhas
7179 the Poisson distribution with mean |ε|≤m|β1|4α+↓|4|≤m|β2.
7185 |π(b) |εN|β1|4α_↓|4N|β2 |πhas the Poisson distribution
7191 with mean |ε|≤m|β1|4α_↓|4|≤m|β2.|'{A3}|≡2|≡7|≡.|9|4[|ε|*/|↔P|
7194 ↔P|\] (|πJ. H. Ahrens.) On most binary computers
7202 there is an e∃cient way to count the number of
7212 1's in a binary word (cf. Section 7.1). Hence
7221 there is a nice way to obtain the binomial distribution
7231 (|εt,|4p) |πwhen |εp|4α=↓|4|f1|d32|), |πsimply
7235 by generating |εt |πrandom bits and counting
7242 the number of |ε|*/|↔O|\|π's.|'!|9|7|πDesign an
7248 algorithm which produces the binomial distribution
7254 (|εt,|4p) |πfor |εarbitrary p, |πusing only a
7261 subroutine for the special case |εp|4α=↓|4|f1|d32|)
7267 |πas a source of random data. [|εHint: |πSimulate
7275 a process which _rst looks at the most signi_cant
7284 bits of |εt |πuniform deviates, then at the second
7293 bit of those deviates whose leading bit is not
7302 su∃cient to determine whether or not their value
7310 is |→|¬W|εp, |πetc.]|'{A3}|≡2|≡8|≡.|9|4[|ε|*/|↔P|↔C|\]
7314 (|πA. J. Walker.) The _nite distribution in exercise
7322 20 could also be generated by the following method:
7331 ``Generate a uniform deviate |εU, |πand set |εN|5|¬L|5|"l6U|
7338 "L. |πIf |εU|4|¬W|4P|βN |πset |εX|4|¬L|4x|βN|βα+↓|β1,
7343 |πotherwise set |εX|4|¬L|4Q|βN, |πwhere (|εP|β0,|4.|4.|4.|4,
7347 |4P|β5)|4α=↓|4(|f256|d31536|), |f499|d31536|),
7349 |f745|d31536|), |f798|d31536|), |f1120|d31536|),
7352 |f1535|d31536|)) |πand (|εQ|β1,|4.|4.|4.|4,|4Q|β5)|4α=↓|4(x|
7354 β1,|4x|β6,|4x|β6,|4x|β3,|4x|β1).'' |πShow that
7357 the general _nite distribution (2) can always
7364 be handled in a similar way, with appropriate
7372 tables (|εP|β0,|4.|4.|4.|4,|4P|βk|βα_↓|β1), (Q|β0,|4.|4.|4.|
7374 4,|4Q|βk|βα_↓|β1).|'{A3}|≡2|≡9|≡.|9|4[|ε|*/HM|↔L|↔C|\]
7376 (|πR. P. Brent.) Develop a method to generate
7384 a random point on the surface of the ellipsoid
7393 de_ned by |¬K|4a|βkx|ur2|)k|)|4α=↓|41, |πwhere
7397 |ε|εa|β1|4|¬R|4|¬O|4|¬O|4|¬O|4|¬R|4a|βn|4|¬Q|40.|'
7398 |H{A6}{H10L12M29}|πThe (|εt|4α+↓|41)|πst record
7401 should be selected with probability (|εn|4α_↓|4m)/(N|4α_↓|4t
7406 ), |πif |εm |πitems have already been selected.
7414 This is the appropriate probability, since of
7421 all the possible ways to choose |εn |πthings
7429 from |εN |πsuch that |εm |πvalues occur in the
7438 _rst |εt, |πexactly|'{A9}|ε|↔a|(N|4α_↓|4t|4α_↓|41|d5n|4α_↓|4
7441 m|4α_↓|41|)|↔s|↔h|↔a|(N|4α_↓|4t|d5n|4α_↓|4m|)|↔s|4α=↓|4|(n|4
7441 α_↓|4m|d2N|4α_↓|4t|)|E|;(1)|?{A9}|πof these select
7446 the (|εt|4α+↓|41)|πst element.|'!|9|7The idea
7451 developed in the preceding paragraph leads immediately
7458 to the following algorithm:|'|Hβ{U0}{H10L12M29}W58320#Comput
folio 182 galley 11
7462 er Programming|9(Knuth/Addison-Wesley)|9f. 182|9ch.
7465 3|9g. 11c|'{A6}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡S
7469 (|εSelection sampling technique).|9|4|πTo select
7473 |εn |πrecords at random from a set of |εN, |πwhere
7483 |ε0|4|¬W|4n|4|¬E|4N.|'{A3}{I1.7H}|≡S|≡1|≡.|9|π[Initialize.]
7485 Set |εt|4|¬L|40, m|4|¬L|40. (|πDuring this algorithm,
7491 |εm |πrepresents the number of records selected
7498 so far, and |εt |πis the total number of input
7508 records we have dealt with.)|'{A3}|≡S|≡2|≡.|9[|πGenerate
7514 |εU.] |πGenerate a random number |εU, |πuniformly
7521 distributed between zero and one.|'{A3}|≡S|≡3|≡.|9[Test.]
7527 If (|εN|4α_↓|4t)U|4|¬R|4n|4α_↓|4m, |πgo to step
7532 S5.|'{A3}|≡S|≡4|≡.|9[Select.] Select the next
7537 record for the sample, and increase |εm |πand
7545 |εt |πby 1. If |εm|4|¬W|4n, |πgo to step S2;
7554 otherwise the sample is complete and the algorithm
7562 terminates.|'{A3}|≡S|≡5|≡.|9[|πSkip.] Skip the
7566 next record (do not include it in the sample),
7575 increase |εt |πby 1, and go to step S2.|'{A12}{IC}|π!|9|7Thi
7584 s algorithm, at _rst glance, may appear to be
7593 unreliable and, in fact, incorrect; but a careful
7601 analysis (see the exercises below) shows that
7608 it is completely trustworthy. It is not di∃cult
7616 to verify that|'{A3}a)|5|5At most |εN |πrecords
7623 are input (we never run o= the end of the _le
7634 before choosing |εn |πitems).|'{A3}b)|9The sample
7640 is completely unbiased; in particular, the probability
7647 that any given element is selected, e.g., the
7655 last element of the _le, is |εn/N.|'{A3}!|9|7|πStatement
7663 (b) is true in spite of the fact that we are
7674 |εnot |πselecting the (|εt|4α+↓|41)|πst item
7679 with probability |εn/N, |πwe select it with the
7687 probability in Eq. (1)*3 This has caused some
7695 confusion in the published literature. Can the
7702 reader explain this seeming contradiction?|'!|9|7(|εNote|*/:
7708 |\|πWhen using Algorithm S, care should be taken
7716 to use a di=erent source of random numbers |εU
7725 |πeach time the program is run, to avoid connections
7734 between the samples obtained on di=erent days.
7741 This can be done, for example, by choosing a
7750 di=erent value of |εX|β0 |πfor the linear congruential
7758 method each time; |εX|β0 |πcould be set to the
7767 current date, or to the last |εX |πvalue generated
7776 on the previous run of the program.)|'!|9|7We
7784 will usually not have to pass over all |εN |πrecords;
7794 in fact, since (b) above says the alst record
7803 is selected with probability |εn/N, |πwe will
7810 terminate the algorithm |εbefore |πconsidering
7815 the last record exactly (1|4α_↓|4|εn/N) |πof
7821 the time. The average number of records considered
7829 when |εn|4α=↓|42 |πis about |f2|d33|)|εN, |πand
7835 the general formulas are given in exercises 5
7843 and 6.|'!|9|7A problem arises if we don't know
7852 the value of |εN |πin advance, since the precise
7861 value of |εN |πis crucial in Algorithm S. Suppose
7870 we want to select |εn |πitems at random from
7879 a _le, without knowing exactly how many are present
7888 in that _le. We could _rst go through and count
7898 the records, then take a second pass to select
7907 them; but it is generally better to sample |εm|4|¬R|4n
7916 |πof the original items on the _rst pass, where
7925 |εm |πis much less than |εN, |πso that only |εm
7935 |πitems must be considered on the second pass.
7943 The trick, of course, is to do this in such a
7954 way that the _nal result is a truly random sample
7964 of the original _le.|'!|9|7|πSince we don't know
7972 when the input is going to end, we must keep
7982 track of a random sample of the input records
7991 seen so far, thus always being prepared for the
8000 end. As we read the input we will construct a
8010 ``reservoir'' which contains only those |εm |πrecords
8017 which have appeared among the previous samples.
8024 The _rst |εn |πrecords always go into the reservoir.
8033 When the (|εt|4α+↓|41)|πst record is being input,
8040 for |εt|4|¬E|4n, |πwe will have in memory a table
8049 of |εn |πindices pointing to those records in
8057 the reservoir which belong to the random sample
8065 we have chosen from the _rst |εt |πrecords. The
8074 problem is to maintain this situation with |εt
8082 |πincreased by one, namely to _nd a new random
8091 sample from among the |εt|4α+↓|41 |πrecords now
8098 known to be present. It is not hard to see that
8109 we should include the new record in the new sample
8119 with probability |εn/(t|4α+↓|41), |πand in such
8125 a case it should replace a random element of
8134 the previous sample.|'!|9|7Thus, the following
8140 procedure does the job:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
8145 |≡R (|εReservoir sampling).|9|4|πTo select |εn
8150 |πrecords at random from a _le of unknown size
8159 |→|¬R|εn, |πgiven |εn|4|¬Q|40. |πAn auxiliary
8164 _le called the ``reservoir'' contains all records
8171 which are condidates for the _nal sample. The
8179 algorithm uses a table of distinct indices |εI[j],
8187 1|4|¬E|4j|4|¬E|4n, |πeach of which points to
8193 one of the records in the reservoir.|'{A3}{I1.8H}|π|≡R|≡1|≡.
8200 |9[Initialize.] Input the _rst |εn |πrecords
8206 and copy them to the reservoir. Set |εI[j]|4|¬L|4j
8214 |πfor |ε1|4|¬E|4j|4|¬E|4n, |πand set |εt|5|¬L|5m|5|¬L|5n.
8219 |πDuring this algorithm, |¬T|εI[1],|4.|4.|4.|4,|4I[n]|¬S
8223 |πpoint to the records in the current sample,
8231 |εm |πis the size of the reservoir, and |εt |πis
8241 the number of input records dealt with so far.)|'
8250 {A3}|≡R|≡2|≡.|9[End of _le?] If there are no
8257 more records to be input, go to step R6.|'{A6}|≡R|≡3|≡.|9[Ge
8266 nerate and test.] Increase |εt |πby 1, then generate
8275 a random integer |εM |πbetween 1 and |εt (|πinclusive).
8284 If |εM|4|¬Q|4n, |πgo to R5.|'{A3}|≡R|≡4|≡.|9[Add
8290 to reservoir.] Copy the next record of the input
8299 _le to the reservoir, increase |εm |πby 1, and
8308 set |εI[M]|5|¬L|5m. (|πThe record previously
8313 pointed to by |εI[M] |πis being replaced in the
8322 sample by the new record.) Go back to R2.|'{A3}|≡R|≡5|≡.|9[S
8331 kip.] Skip over the next record of the input
8340 _le (do not include it in the reservoir), and
8349 return to step R2.|'{A3}|≡R|≡6|≡.|9[Second pass.]
8355 Sort the |εI |πtable entries so that |εI[1]|4|¬W|4|¬O|4|¬O|4
8362 |¬O|4I[n]; |πthen go through the reservoir, copying
8369 the records with these indices into the output
8377 _le which is to hold the _nal sample.|'{A12}{IC}!|9|7Algorit
8385 hm R is due to Alan G. Waterman. The reader may
8396 wish to work out the example of its operation
8405 which appears in exercise 9.|'!|9|7If the records
8413 are su∃ciently short, it is of course unnecessary
8421 to have a reservoir at all; we can keep the |εn
8432 |πrecords of the current sample in memory at
8440 all times (see exercise 10).|'!|9|7|πThe natural
8447 question to ask about Algorithm R is, ``What
8455 is the expected size, of the reservoir?'' Exercise
8463 11 shows that the average value of |εm |πis exactly
8473 |εn(1|4α+↓|4H|βN|4α_↓|4H|βn); |πthis is approximately
8477 |εn{H12}({H10}1|4α+↓|4|πln(|εN/n){H12}){H10}.
8478 So if |εN/n|4α=↓|41000, |πthe reservoir will
8484 contain only a{U0}{H10L12M29}W58320#Computer
folio 186 galley 12
8487 Programming|9(Knuth/Addison-Wesley)|9f. 186|9ch.
8489 3!g. 12c|'{A6}!|9|7|πNote that Algorithms S and
8496 R can be used to obtain samples for several independent
8506 categories simultaneously. For example, if we
8512 have a large _le of names and addresses of U.S.
8522 residents, we could pick random samples of exactly
8530 10 people from each of the 50 states without
8539 making 50 passes through the _le, and without
8547 _rst sorting the _le by state.|'!|9|7Algorithm
8554 S and a number of other sampling techniques are
8563 discussed in a paper by C. T. Fan, Mervin E.
8573 Muller, and Ivan Rezucha, |εJournal of the American
8581 Statistical Association |π|≡5|≡7 (1962), 387<402.
8586 The method was independently discovered by T.
8593 G. Jones, |εCACM |π|≡5 (1962), 343.|'!|9|7The
8600 |εsampling problem, |πas described here, can
8606 be regarded as the computation of a random |εcombination,
8615 |πaccording to the conventional de_nition of
8621 combinations of |εN |πthings taken |εn |πat a
8629 time (see Section 1.2.6). Now let us consider
8637 the problem of computing a random |εpermutation
8644 |πof |εt |πobjects; we will call this the |εshu+ing
8653 problem, |πsince shu+ing a deck of cards is nothing
8662 more than subjecting it to a random permutation.|'
8670 !|9|7A moment's re⊗ection is enough to convince
8677 oneself that the approaches people traditionally
8683 use to shu+e cards are miserably inadequate;
8690 there is no hope of obtaining each of the |εt*3
8700 |πpermutations with anywhere near equal probability
8706 by such methods. It has been said that expert
8715 bridge players make use of this fact when deciding
8724 whether or not to ``_nesse.''|'!|9|7|πWhen |εt
8731 |πis small, we can obtain a random permutation
8739 very quickly by generating a random integer between
8747 1 and |εt*3. |πFor example, when |εt|4α=↓|44,
8754 |πa random number between 1 and 24 su∃ces to
8763 select a random permutation from a table of all
8772 possibilities. But for alrge |εt, |πit is necessary
8780 to be more careful if we want to claim that each
8791 permutation is equally likely, since |εt*3 |πis
8798 much larger than the accuracy of individual random
8806 numbers.|'!|9|7A suitable shu+ing procedure can
8812 be obtained by recalling Algorithm 3.3.2P, which
8819 gives a simple correspondence between each of
8826 the |εt*3 |πpossible permutations and a sequence
8833 of numbers (|εc|β1,|4c|β2,|4.|4.|4.|4,|4c|βt|βα_↓|β1)
8836 |πwith |ε0|4|¬E|4c|βj|4|¬E|4j. |πIt is easy to
8842 compute such a set of numbers at random, and
8851 we can use the correspondence to produce a random
8860 permutation.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
8862 |≡P (|εShu∪ing).|9|4|πLet |εX|β1, X|β2, |4.|4.|4.|4,|4X|βt
8867 |πbe a set of |εt |πnumbers to be shu+ed.|'{A3}{I1.8H}|≡P|≡1
8876 |≡.|9[Initialize.] Set |εj|5|¬L|5t.|'{A3}|π|≡P|≡2|≡.|9[Gener
8879 ate |εU.] |πGenerate a random number U, |πuniformly
8887 distributed between zero and one.|'{A3}|≡P|≡3|≡.|9[Exchange.
8892 ] Set |εk|5|¬L|5|"ljU|"L|4α+↓|41. (|πNow |εk
8897 |πis a random integer between 1 and |εj). |πExchange
8906 |εX|βk|5|≠m|5X|βj.|'{A3}|π|≡P|≡4|≡.|9[Decrease
8908 |εj.] |πDecrease |εj |πby 1. If |εj|4|¬Q|41,
8915 |πreturn to step P2.|'{A12}{IC}!|9|7This algorithm
8921 was _rst published by L. E. Moses and R. V. Oakford,
8932 in |εTables of Random Permutations (|πStanford
8938 University Press, 1963), and by R. Durstenfeld,
8945 |εCACM |π|≡7 (1964), 420. It can be used also
8954 to obtain random combinations (see exercise 14).|'
8961 {A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'{A12}{H9L11}|5|5|≡1|≡.|9|4
8962 [|ε|*/M|↔O|↔P|\] |πExplain Eq. (1).|'{A3}|5|5|≡2|≡.|9|4[|ε|*/|
8966 ↔P|↔c|\] |πProve that Algorithm S never tries
8973 to read more than |εN |πrecords of its input
8982 _le.|'{A3}|5|5|≡3|≡.|9|4[|ε|*/|↔P|↔P|π|\] The
8985 (|εt|4α+↓|41)|πst item in Algorithm S is selected
8992 with probability (|εn|4α_↓|4m)/(N|4α_↓|4t), not
8996 n/N, |πyet the text claims that the sample is
9005 unbiased so each item should be selected with
9013 the |εsame |πprobability*3 How can both of these
9021 statements be true?|'{A3}|5|5|≡4|≡.|9|4[|ε|*/M|↔P|↔L|\]
9025 |πLet |εp(m,|4t) |πbe the probability that exactly
9032 |εm |πitems are selected from among the _rst
9040 |εt |πin the selection sampling technique. Show
9047 directly from Algorithm S that|'{A9}|εp(m,|4t)|4α=↓|4|↔a|(t|
9052 d5m|)|↔s|↔a|(N|4α_↓|4t|d5n|4α_↓|4m|)|↔s|↔h|↔a|(N|d5n|)|↔s,!!
9052 |πfor!!|ε0|4|¬E|4t|4|¬E|4N.|;{A9}|5|5|≡5|≡.|9|4[|ε|*/M|↔P|↔M|
9053 \] |πWhat is the average value of |εt |πwhen
9062 Algorithm S terminates? (In other words, how
9069 many of the |εN |πrecords have passed before
9077 the sample is complete?)|'{A3}|5|5|≡6|≡.|9|4[|ε|*/M|↔p|↔M|\]
9082 |πWhat is the standard deviation of the value
9090 computed in the previous exercise?|'{A3}|5|5|≡7|≡.|9|4[|ε|*/M
9095 |↔P|↔C|\] |πProve that any |εgiven |πchoice of
9102 |εn |πrecords from the set of |εN |πis obtained
9111 by Algorithm S with probability 1/(|ε|fN|d5n|)).
9117 |πTherefore the sample is completely unbiased.|'
9123 |5|5|≡8|≡.|9|4[|ε|*/M|↔M|↔o|\] |πAlgorithm S computes
9127 one uniform deviate for each input record it
9135 handles. Find a more e∃cient way to determine
9143 the number of input records to skip before the
9152 _rst is selected, assuming that |εN/n |πis rather
9160 large. (We could iterate this process to select
9168 the remaining |εn|4α_↓|41 |πrecords.)|'|5|5|≡9|≡.|9|4[|ε|*/|↔
9172 O|↔P|\] |πLet |εn|4α=↓|43. |πIf Algorithm R is
9179 applied to a _le containing 20 records numbered
9187 1 thru 20, and if the random numbers generated
9196 in step R3 are respectively|'{A6}4,|41,|46,|47,|45,|43,|45,|
9201 411,|411,|43,|47,|49,|43,|411,|44,|45,|44,|;{A6}which
9203 records go into the reservoir? Which are in the
9212 _nal sample?|'{A3}|≡1|≡0|≡.|9|4[|ε|*/|↔O|↔C|\]
9215 |πModify Algorithm R so that the reservoir is
9223 eliminated, assuming that the |εn |πrecords of
9230 the current sample can be held in memory.|'{A3}|≡1|≡1|≡.|9|4
9238 [|ε|*/M|↔P|↔C|\] |πLet |εp|βm |πbe the probability
9244 that exactly |εm |πelements are put into the
9252 reservoir during the _rst pass of Algorithm R.
9260 Determine the generating function |εG(z)|4α=↓|4|¬K|βmp|βmz|g
9264 m, |πand _nd the mean and standard deviation.
9272 (Use the ideas of Section 1.2.10.)|'{A3}|≡1|≡2|≡.|9|4[|ε|*/M|
9278 ↔P|↔o|\] |πThe gist of Algorithm P is that any
9287 permutation |ε|≤p |πcan be uniquely written as
9294 a product of transpositions in the form |ε|≤p|4α=↓|4(a|βtt)|
9301 4.|4.|4.|4(a|β33)(a|β22), |πwhere |ε1|4|¬E|4a|βj|4|¬E|4j
9304 |πfor |εt|4|¬R|4j|4|¬Q|41. |πProve that there
9309 is also a unique representation of the form |ε|≤p|4α=↓|4(b|β
9317 22)(b|β33)|4.|4.|4.|4(b|βtt), |πwhere |ε1|4|¬E|4b|βj|4|¬E|4j
9319 |πfor |ε1|4|¬W|4j|4|¬E|4t, |πand give an algorithm
9326 which computes the |εb'|πs from the |εa'|πs in
9334 |εO(t) |πst{U0}{H10L12M29}W58320#Computer Programming|9(Knut
folio 190 galley 13
9336 h/Addison-Wesley)|9f. 190|9ch. 3|9g. 13c|'{A6}{H9L11}|≡1|≡3|
9340 ≡.|9|4[|ε|*/M|↔P|↔L|\] |π(S. W. Golomb.) A common
9346 procedure for card shu+ing is to divide the deck
9355 into two parts as equal as possible, and to ``ri+e''
9365 them together. (See the discussion of card-playing
9372 etiquette in Hoyle's rules of card games; we
9380 read, ``A shu+e of this sort should be made about
9390 three times to mix the cards thoroughly.'') Consider
9398 a deck of |ε2n|4α_↓|41 |πcards |εX|β1,|4X|β2,|4.|4.|4.|4,|4X
9403 |β2|βn|βα_↓|β1; |πa ``perfect shu+e'' |εs |πdivides
9409 this deck into |εX|β1,|4X|β2,|4.|4.|4.|4,|4X|βn
9413 |πand |εX|βn|βα+↓|β1,|4.|4.|4.|4,|4X|β2|βn|βα_↓|β1
9415 |πand perfectly interleaves them, to obtain |εX|β1,
9422 X|βn|βα+↓|β1, X|β2, X|βn|βα+↓|β2,|4.|4.|4.|4,
9425 X|β2|βn|βα_↓|β1,|4X|βn. |πThe ``cut'' operation
9429 |εc|gj |πchanges |εX|β1, X|β2,|4.|4.|4.|4, X|β2|βn|βα_↓|β1
9434 |πinto |εX|βj|βα+↓|β1,|4.|4.|4.|4, X|β2|βn|βα_↓|β1,|4X|β1,|4
9436 .|4.|4.|4,|4X|βj. |πShow that by combining perfect
9442 shu+es and cuts, at most (|ε2n|4α_↓|41)(2n|4α_↓|42)
9448 |πdi=erent arrangements of the deck are possible,
9455 if |εn|4|¬Q|41.|'|≡1|≡4|≡.|9|4[|ε|*/|↔L|↔c|\]
9458 |π(Ole-Johan Dahl.) If |εX|βj|4α=↓|4j |πat the
9464 start of Algorithm P, and if we terminate that
9473 algorithm when |εj |πreaches the value |εt|4α_↓|4n,
9480 |πthe sequence |εX|βt|βα_↓|βn|βα+↓|β1,|4.|4.|4.|4,|4X|βt
9483 |πis a random permutation of a random combination
9491 of |εn |πelements. Show how to simulate the e=ect
9500 of this procedure using only |εO(n) |πcells of
9508 memory.|'{A24}{H10L12}|≤⊂|∨3|∨.|∨5|∨.|9|∨W|∨H|∨A|∨T|9|∨I|∨S|
9509 9|∨A|9|∨R|∨A|∨N|∨D|∨O|∨M|9|∨S|∨E|∨Q|∨U|∨E|∨N|∨C|∨E|∨*3|'
9510 {A12}|≡A|≡.|9|≡I|≡n|≡t|≡r|≡o|≡d|≡u|≡c|≡t|≡o|≡r|≡y
9511 |∨r|∨e|∨m|∨a|∨r|∨k|∨s|∨.|9|4We have seen in this
9516 chapter how to generate sequences|'{A9}|ε|¬DU|βn|¬F|4α=↓|4U|
9521 β0,|4U|β1,|4U|β2,|4.|4.|4.|E|;(1)|?{A9}|πof real
9525 numbers in the range |ε0|4|¬E|4U|βn|4|¬W|41,
9530 |πand we have called them ``random'' sequences
9537 even though they are completely deterministic
9543 in character. To justify this terminology, we
9550 claimed that the numbers ``behave as if they
9558 are truly random.'' Such a statement may be satisfactory
9567 for practical purposes (at the present time),
9574 but it sidesteps a very important philosophical
9581 and theoretical question: Precisely what do we
9588 mean by ``random behavior''? A quantitative de_nition
9595 is needed. It is undesirable to talk about concepts
9604 that we do not really understand, especially
9611 since many apparently paradoxical statements
9616 can be made about random numbers.|'!|9|7The mathematical
9624 theory of probability and statistics carefully
9630 avoids answering the question; it refrains from
9637 making absolute statements, and instead expresses
9643 everything in terms of how much |εprobability
9650 |πis to be attached to statements involving random
9658 sequences of events. The axioms of probability
9665 theory are set up so that abstract probabilities
9673 can be computed readily, but nothing is said
9681 about what probability really signi_es, or how
9688 this concept can be applied meaningfully to the
9696 actual world. In the book |εProbability, Statistics,
9703 and Truth (|πMacmillan, 1957), R. von Mises discusses
9711 this situation in detail, and presents the view
9719 that a proper de_nition of probability depends
9726 on obtaining a proper de_nition of a rafdbmm
9734 s*?*?*?oper de_nition of a random sequence.|'!|9|7Let
9741 us paraphrase here some statements made by two
9749 of the many authors who have commented on the
9758 subject.|'{A12}{I1.7l}|εD. H. Lehmer (|*/|↔O|↔m|↔C|↔O|\)|*/:
9763 |\|π``A random sequence is a vague notion embodying
9771 the idea of a sequence in which each term is
9781 unpredictable to the uninitiated and whose digits
9788 pass a certain number of tests, traditional with
9796 statisticians and depending somewhat on the uses
9803 to which the sequence is to be put.``|'{A6}|εJ.
9812 N. Franklin (|*/|↔O|↔m|↔o|↔P|\)|*/: |\|π``The sequence
9817 (1) is random if it has every property that is
9827 shared by all in_nite sequences of independent
9834 samples of random variables from the uniform
9841 distribution.''|'{IC}{A12}!|9|7|πFranklin's statement
9844 essentially generalizes Lehmer's to say that
9850 the sequence must satisfy |εall |πstatistical
9856 tests. His de_nition is not completely precise,
9863 and we will see later that a reasonable interpretation
9872 of his statement leads us to conclude that there
9881 is no such thing as a random sequence*3 Let us
9891 begin with Lehmer's less restrictive statement
9897 and attempt to make |εit |πprecise. What we really
9906 want is a relatively short list of mathematical
9914 properties, each of which is satis_ed by our
9922 intuitive notion of a random sequence; furthermore,
9929 this list is to be complete enough so that we
9939 are willing to agree that |εany |πsequence satisfying
9947 these properties is ``random.'' In this section,
9954 we will develop what seems to be an adequate
9963 de_nition of randomness according to these criteria,
9970 although many interesting questions remain to
9976 be answered.|'!|9|7Let |εu |πand |εv |πbe real
9984 numbers, |ε0|4|¬E|4u|4|¬W|4v|4|¬E|41. |πIf |εU
9988 |πis a random variable which is uniformly distributed
9996 between 0 and 1, the probability that |εu|4|¬E|4U|4|¬W|4v
10004 |πis equal to |εv|4α_↓|4u. |πFor example, the
10011 probability that |f1|d33|)|4|¬E|4U|4|¬W|4|f2|d33|)
10014 is |f1|d33|). |πHow can we translate this property
10022 of the single number |εU |πinto a property of
10031 the in_nite sequence |εU|β0,|4U|β1,|4U|β2,|4.|4.|4.|4?
10035 |πThe obvious answer is to count how many times
10044 |εU|βn |πlies between |εu |πand |εv, |πand the
10052 average number of times should equal |εv|4α_↓|4u.
10059 |πOur intuitive idea of probability is based
10066 in this way on the frequency of occurrence. More
10075 precisely, let |εv(n) |πbe the number of values
10083 of |εj, 0|4|¬E|4j|4|¬W|4n, |πsuch that |εu|4|¬E|4U|βj|4|¬W|4
10088 v; |πwe want the ratio |ε|≤n(n)/n |πto approach
10096 the value |εv|4α_↓|4u |πas |εn |πapproaches in_nity:|'
10103 {A9}|uwlim|uc|)|.|εn|¬M|¬X|)|4|≤n(n)/n|4α=↓|4v|4α_↓|4u.|;
10104 {A9}|πIf this condition holds for all choices
10111 of |εu |πand |εv, |πthe sequence is said to be
10121 |εequidistributed.|'!|9|7|πLet |εS(n) |πbe a
10126 statement about the integer |εn |πand the sequence
10134 |εU|β1,|4U|β2,|4.|4.|4.|4; |πfor example, |εS(n)
10138 |πmight be the statement considered above, namely
10145 ``u|4|¬E|4U|βn|4|¬W|4v.'' |πWe can generalize
10150 the idea used in the preceding paragraph to de_ne
10159 ``the probability that |εS(n) |πis true'' with
10166 respect to a particular in_nite sequence: Let
10173 |ε|≤n(n) |πbe the number of values of |εj, 0|4|¬E|4j|4|¬W|4n
10181 , |πsuch that |εS(j) |πis true.|'{A12}|≡D|≡e|≡_|≡n|≡i|≡t|≡i|
10187 ≡o|≡n |≡A|≡.|9|4|εWe say |πPr{H12}({H10}|εS(n){H12}){H10}|4α
10190 =↓|4|≤l, if |πlim|ε|βn|β|¬M|β|¬X|4|≤n(n)/n|4α=↓|4|≤l.
10193 (|πRead, ``The probability that |εS(n) |πis true
10200 is |ε|≤l, |πif the limit as |εn |πtends to in_nity
10210 of |ε|≤n(n)/n |πis |ε|≤l.'')|'{A12}|πIn terms
10216 of this notation, the sequence |εU|β0,|4U|β1,|4.|4.|4.
10222 |πis equidistributed if and only if Pr(|εu|4|¬E|4U|βn|4|¬W|4
10228 v)|4α=↓|4v|4α_↓|4u, |πfor all real numbers |εu,
10234 v |πwith |ε0|4|¬E|4u|4|¬W|4v|4|¬E|41.|'!|9|7|πA
10238 sequence may be equidistributed without being
10244 random. For example, if |εU|β0, U|β1,|4.|4.|4.
10250 |πand |εV|β0,|4V|β1,|4.|4.|4. |πare equidistributed
10254 sequences, it is not hard to show that the sequence|'
10264 {A9}|εW|β0,|4W|β1,|4W|β2,|4W|β3,|4.|4.|4.|4α=↓|4|f1|d32|)U|β
10264 0,|4|f1|d32|)|4α+↓|4|f1|d32|)V|β0,|4|f1|d32|)U|β1,|4|f1|d32|
10264 )|4α+↓|4|f1|d32|)V|β1,|4.|4.|4.|E|;(3)|?{A9}|πis
10267 also equidistributed, since the sequence |ε|f1|d32|)U|βo,|4|
10272 f1|d32|)U|β1,|4.|4.|4. |πis equidistributed between
10276 0 and |f1|d32|), |πwhile the alternate terms,
10283 |f1|d32|)|4α+↓|4|f1|d32|)|εV|β0,|4|f1|d32|)|4α+↓|4|f1|d32|)V
10283 |β1,|4.|4.|4.|4, |πare equidistributed between
10287 |f1|d32|) and 1. In the sequence of |εW|1'|πs,
10295 a value less than |f1|d32|) is always followed
10303 by a value greater than or equal to |f1|d32|),
10312 and conversely; hence that sequence is not random
10320 by any reasonable de_nition. A stronger property
10327 than equidistribution is needed.|'!|9|7A natural
10333 generalization of the equidistribution property,
10338 which removes the objection stated in the preceding
10346 paragraph, is to consider adjacent pairs of numbers
10354 of our sequence. We can require the sequence
10362 to satisfy the condition|'{A9}Pr(|εu|β1|4|¬E|4U|βn|4|¬W|4v|β
10366 1!!|πand!!|εu|β2|4|¬E|4U|βn|βα+↓|β1|4|¬W|4v|β2)|4α=↓|4(v|β1|
10366 4α_↓|4u|β1)(v|β2|4α_↓|4u|β2)|E|;(4)|?{A9}|πfor
10369 any four numbers |εu|β1, v|β1, u|β2, v|β2 |πwith
10377 |ε0|4|¬E|4u|β1|4|¬W|4v|β1|4|¬E|41,|40|4|¬E|4u|β2|4|¬W|4v|β2|
10377 4|¬E|41. |πIn general, for any positive integer
10384 |εk |πwe can require our sequence to be |εk-distributed
10393 |πin the following sense:|'{A12}|≡D|≡e|≡_|≡n|≡i|≡t|≡i|≡o|≡n
10398 |≡B|≡.|9|4|εThe sequence (|*/|↔O|\) is said to
10404 be k-distributed if|'{A9}|πPr(|εu|β1|4|¬E|4U|βn|4|¬W|4v|β1,|
10407 4.|4.|4.|4,|4u|βk|4|¬E{U0}{H10L12M29}W58320#Computer
folio 194 galley 14
10408 Programming|9(Knuth/Addison-Wesley)|9f. 194|9ch.
10410 3|9g. 14c|'{A6}!|9|7|πAn equidistributed sequence
10415 is a 1-distributed sequence. Note that if |εk|4|¬Q|4,
10423 |πa |εk-|πdistributed sequence is always (|εk|4α_↓|41)-|πdis
10428 tributed, since we may set |εu|βk|4α=↓|40 |πand
10435 |εv|βk|4α=↓|41 |πin Eq. (5). Thus, in particular,
10442 any sequence which is known to be 4-distributed
10450 must also be 3-distributed, 2-distributed, and
10456 equidistributed. We can investigate the largest
10462 |εk |πfor which a given sequence is |εk-|πdistributed;
10470 and this leads us to formulate|'{A12}|≡D|≡e|≡_|≡n|≡i|≡t|≡i|≡
10476 o|≡n |≡C|≡.|9|4|εA sequence is said to be |¬X-distributed
10484 if it is k-distributed for all positive integers
10492 k.|'{A12}!|9|7|πSo far we have considered ``[0,|41)
10499 sequences,'' i.e., sequences of real numbers
10505 lying between zero and one. The same ideas apply
10514 to integer-valued sequences; let us say a sequence
10522 |ε|¬DX|βn|¬F|4α=↓|4X|β0,|4X|β1,|4X|β2,|4.|4.|4.
10523 |πis a ``|εb-|πary sequence'' if each |εX|βn
10530 |πis one of the integers 0,|41,|4.|4.|4.|4,|4b|4α_↓|41.
10536 |πThus, a 2-ary (binary) sequence is a sequence
10544 of zeros and ones.|'!|9|7A |εk-|πdigit ``|εb-|πary
10551 number'' |εx|β1x|β2|4.|4.|4.|4x|βk |πis an ordered
10556 set of |εk |πintegers, where |ε0|4|¬E|4x|βj|4|¬W|4b
10562 |πfor |ε1|4|¬E|4j|4|¬E|4k.|'{A12}|π|≡D|≡e|≡_|≡n|≡i|≡t|≡i|≡o|
10564 ≡n |≡D|≡.|9|4|εA b-ary sequence is said to be
10572 k-distributed if|'{A9}|πPr(|εX|βnX|βn|βα+↓|β1|4.|4.|4.|4X|βn
10574 |βα+↓|βk|βα_↓|β1|4α=↓|4x|β1x|β2|4.|4.|4.|4x|βk)|4α=↓|41/b|gk
10574 |E|;(6)|?{A9}for all b-ary numbers x|β1x|β2|4.|4.|4.|4x|βk.|
10580 '{A12}!|9|7|πIt is clear from this de_nition
10587 that if |εU|β0,|4U|β1,|4.|4.|4. |πis a |εk-|πdistributed
10593 [0,|41) sequence, then the sequence |↔l|εbU|β0|↔L,
10599 |↔lbU|β1|↔L,|4.|4.|4. |πis a |εk-|πdistributed
10603 |εb-|πary sequence. (For if we set |εu|βj|4α=↓|4x|βj/b,
10610 v|βj|4α=↓|4(x|βj|4α+↓|41)/b, X|βn|4α=↓|4|↔lbU|βn|↔L,
10612 |πEq. (5) becomes Eq. (6).{H12}){H10} Furthermore,
10618 every |εk-|πdistributed |εb-|πary sequence, for
10623 |εk|4|¬Q|41, |πis also (|εk|4α_↓|41)-|πdistributed:
10627 we add together the probabilities for the |εb-|πary
10635 numbers |εx|β1|4.|4.|4.|4x|βk|βα_↓|β10, x|β1|4.|4.|4.|4x|βk|
10637 βα_↓|β11,|4.|4.|4.|4, x|β1|4.|4.|4.|4x|βk|βα_↓|β1(b|4α_↓|41)
10638 |πto obtain|'{A9}|πPr(|εX|βn|4.|4.|4.|4X|βn|βα+↓|βk|βα_↓|β2
10641 |4α=↓|4x|β1|4.|4.|4.|4x|βk|βα_↓|β1)|4α=↓|41/b|gk|gα_↓|g1.|;
10642 {A9}(|πProbabilities for disjoint events are
10647 additive; see exercise 5.) It therefore is natural
10655 to speak of an |¬X-distributed |εb-|πary sequence,
10662 as in De_nition C above.|'!|9|7|πThe representation
10669 of a positive real number in the radix-|εb |πnumber
10678 system may be regarded as a |εb-|πary sequence;
10686 for example, |ε|≤p |πcorresponds to the 10-ary
10693 sequence 3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8,
10706 9,|4.|4.|4.|4. |πIt has been conjectured that
10712 this sequence is |¬X-distributed, but nobody
10718 has yet been able to prove that it is even 1-distributed.|'
10729 !|9|7|πLet us analyze these concepts a little
10736 more closely in the case when |εk |πequals a
10745 million. A binary sequence which is 1000000-distributed
10752 is going to have runs of a million zeros in a
10763 row*3 Similarly, a [0,|41) sequence which is 1000000-distribu
10770 ted is going to have runs of a million consecutive
10780 values each of which is less than |f1|d32|).
10788 It is true that this will happen only (|f1|d32|))|g1|g0|g0|g
10796 0|g0|g0|g0 of the time, on the average, but the
10805 fact is that it |εdoes |πhappen. Indeed, this
10813 phenomenon will occur in any truly random sequence,
10821 using our intuitive notion of ``truly random.''
10828 One can easily imagine that such a situation
10836 will have a drastic e=ect if this set of a million
10847 ``truly random'' numbers is being used in a computer-simulat
10855 ion experiment; there would be good reason to
10863 complain about the random-number generator*3 However,
10869 if we have a sequence of numbers which never
10878 has runs of a million consecutive |εU'|πs less
10886 than |f1|d32|), the sequence is not random, and
10894 it will not be a suitable source of numbers for
10904 other conceivable applications which use extremely
10910 long blocks of |εU'|πs as input. In summary,
10918 |εa truly random sequence will exhibit local
10925 nonrandomness|*/; |\|πlocal nonrandomness is necessary
10930 in some applications, but it is disastrous in
10938 others. We are forced to conclude that |εno sequence
10947 of ``random'' numbers can be adequate for every
10955 application.|'!|9|7|πIn a similar vein, one may
10962 argue that there is no way to judge whether a
10972 |ε⊂nite |πsequence is random or not; any particular
10980 sequence is just as likely as any other one.
10989 These facts are de_nitely stumbling blocks if
10996 we are ever to have a useful de_nition of randomness,
11006 but they are not really cause for alarm. It is
11016 still possible to give a de_nition for the randomness
11025 of in_nite sequences of real numbers in such
11033 a way that the corresponding theory (viewed properly)
11041 will give us a great deal of insight concerning
11050 the ordinary _nite sequences of rational numbers
11057 which are actually generated on a computer. Furthermore,
11065 we shall see later in this section that there
11074 are several plausible de_nitions of randomness
11080 for _nite sequences are in fact available.|'{A12}|≡B|≡.|9|¬X
11087 |≡-|≡d|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡e|≡d |≡s|≡e|≡q|≡u|≡e|≡n|≡c|≡
11088 e|≡s|≡.|9|4Let us now undertake a brief study
11095 of the theory of |¬X-distributed sequences. To
11102 describe the theory adequately, we will need
11109 to use a little ``higher mathematics,'' so we
11117 assume in the remainder of this subsection that
11125 the reader knows the material ordinarily taught
11132 in an ``advanced calculus'' course.|'|Ha*?*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!
11137